// Billiard words of length <= n in m dimensions, Fred Lunnon, August 2010  
// Language represented by sorted list of integers 
// [No strings ... and no sets please, we're British!] 

bas := 10; // base = 10 or max(m+1) (global preset constant)  

// Kronecker delta 
function krodel (x, y) return x eq y select 1 else 0; end function; 

// Nec and suff test for billiard word via inequalities & linear programming 
function viable (word_int, n, m) 
  local i,j,k,s,t,lhs,rhs,rel,obj,R,cod,word_str; 
  //global bas; 
  RF := RealField(); 
  // s[j] = (integer) symbol in position j of word 
  word_str := IntegerToString(word_int, bas); 
  s := [StringToInteger(word_str[i]) : i in [1..n]]; 
  // t[j] = number previous occurrences of symbol s[j] in word, extended L & R  
  t := [0 : k in [1..m]] cat 
    [&+[krodel(s[j], s[i]) : i in [1..j]] : j in [1..n]] cat  
    [1 + &+[krodel(k, s[i]) : i in [1..n]] : k in [1..m]]; 

  lhs := Matrix(RF, 2*m+n-1, 2*m, [] 
    cat [<k, s[n], 1> : k in [1..m]] 
    cat [<k, s[n]+m, t[n+m]> : k in [1..m]] 
    cat [<k, k, -1> : k in [1..m]] 
    cat [<k, k+m, -t[n+k+m]> : k in [1..m]] 

    cat [<j+m, s[j], 1> : j in [1..n-1]] 
    cat [<j+m, s[j]+m, t[j+m]> : j in [1..n-1]] 
    cat [<j+m, s[j+1], -1> : j in [1..n-1]] 
    cat [<j+m, s[j+1]+m, -t[j+1+m]> : j in [1..n-1]] 

    cat [<k+m+n-1, k, 1> : k in [1..m]] 
    cat [<k+m+n-1, k+m, t[1-k+m]> : k in [1..m]] 
    cat [<k+m+n-1, s[1], -1> : k in [1..m]] 
    cat [<k+m+n-1, s[1]+m, -t[1+m]> : k in [1..m]] 
    ); 

  rhs := Matrix(RF, 2*m+n-1, 1, [-1 : i in [1..2*m+n-1]]); // avoid < relations 
  rel := Matrix(RF, 2*m+n-1, 1, [-1 : i in [1..2*m+n-1]]); // spec <= relations 
  obj := Matrix(RF, 1, 2*m, [0 : i in [1..2*m]]); // dummy objective 

  _,cod := MaximalSolution(lhs, rel, rhs, obj); // 0,2 if solution,infeasible 
  return cod eq 0 select true else cod eq 2 select false else cod; end function; 

// Extend m-symbol n-word with m projections to (m-1)-symbol words: 
//   l = max symbol used in word; code,n,l = zero on entry; set codbas > m. 
//   prolangs contains all (m-1)-symbol words of length <= n. 
//   Inequality test omitted when single possible continuation; 
//   permutations suppressed; maybe check new suffix among n-words? 
// Build option to return languages: 
//   prolang[m] = join of m-symbol lang[n] for n no greater 
procedure extend (oldproj, oldword, maxlen, n, m, l, build, ~totals, ~langs, ~prolang) 
  local i,j,k,s,newproj,newword,kset; 
  //global totals,langs,prolang,bas; 
  totals[n+1] +:= (l eq 0) select 1 else &*[m+1-i : i in [1..l]]; 
    // total + m...(m-l+1) --- & fails if range empty! 
  if build then Append(~langs[n+1], oldword); end if; 
  if n lt maxlen then 
    newproj := [0 : k in [1..m]]; // extended projections 

    // Locate all extensions to symbol k with good prime projections; 
    //   latter already reduced modulo permutation! 
    kset := [k : k in [1..Min(m,l+1)]]; // max symbol increment at most unity  
    for k in kset do 
      for i := 1 to m do if k notin kset then break; end if; 
        // check (m-1)-symbol projections extended by k-1 or k 
        if i lt k then newproj[i] := oldproj[i]*bas + k-1; 
        elif i gt k then newproj[i] := oldproj[i]*bas + k; 
        else newproj[i] := oldproj[i]; end if; 
        if i ne k and newproj[i] notin prolang[m-1+1] then Exclude(~kset, k); 
          // search by bisection slower in practice! 
        end if; end for; 
    end for; 

    // Extend sequence to every remaining symbol k; 
    // solve inequalities only if several branches; assert never none! 
    for k in kset do 
      for i := 1 to m do 
        if i lt k then newproj[i] := oldproj[i]*bas + k-1; 
        elif i gt k then newproj[i] := oldproj[i]*bas + k; 
        else newproj[i] := oldproj[i]; end if; end for; 
      newword := oldword*bas + k; 
      if #kset eq 1 or viable(newword, n+1, m) then // COR here !! 
        extend(newproj, newword, maxlen, n+1, m, Max(k,l), build, 
          ~totals, ~langs, ~prolang); end if; 
    end for; end if; end procedure; 

// Main program: generate lower dimensions m first for projection; 
// depth-first generation of words up to length n 
  maxsym := 4; n := 10; // dimension and length 
  prolang := [[0]] cat [[] : j in [1+1..maxsym+1]]; 
  for m := 1 to maxsym do 
    totals := [0 : j in [0+1..n+1]]; // count of words by length 
    proword := [0 : k in [1+1..m+1]]; // projections of current word 
    langs := [[0]] cat [[] : j in [1+1..n+1]]; 
      // words by length mod perm: initial empty word 
    time extend(proword, 0, n, 0, m, 0, true, ~totals, ~langs, ~prolang); 
    prolang[m+1] := &cat [langs[j] : j in [0+1..n+1]]; // old langs union 
    print m, n, totals; end for; 

/* Output 
Time: 0.000
1 10 [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]
Time: 0.080
2 10 [ 1, 2, 4, 8, 14, 24, 36, 54, 76, 104, 136 ]
Time: 0.460
3 10 [ 1, 3, 9, 27, 75, 189, 447, 951, 1911, 3621, 6513 ]
Time: 2.780
4 10 [ 1, 4, 16, 64, 244, 856, 2776, 8356, 23032, 59200, 142624 ]
*/