// 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 ] */