// 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 in [1..m]] cat [ : k in [1..m]] cat [ : k in [1..m]] cat [ : k in [1..m]] cat [ : j in [1..n-1]] cat [ : j in [1..n-1]] cat [ : j in [1..n-1]] cat [ : j in [1..n-1]] cat [ : k in [1..m]] cat [ : k in [1..m]] cat [ : k in [1..m]] cat [ : 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 := 3; n := 15; // 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; // langs[n][i] holds the i-th word of length n, encoded as a decimal number, // and excluding permutations; bas := 10; // base = 10 or max(m+1) (global preset constant) // Reduce words in langs[n+2] mod 10, convert to multiset; // special[k][n] := number occurring k times. Increment by perms! m := 3; maxlen := 15; baspow := [10^(j-1) : j in [0+1..maxlen+1]]; special := [[0 : j in [0+1..maxlen+1]] : k in [1..m]]; // counts of words by speciality & length for n := 1 to maxlen do i := 0; while i lt #langs[n+1] do i +:= 1; word := Floor(langs[n+1][i] / bas); // count alphabet l := (n eq 1) select 0 else Max([Floor(word / baspow[j]) mod bas : j in [0+1..n-1]]); // count speciality k of prefix k := Max(m-l-1,0); while i le #langs[n+1] and Floor(langs[n+1][i] / bas) eq word do i +:= 1; k +:= 1; end while; // count isomorphs special[k][n] +:= (l eq 0) select 1 else &*[m+1-j : j in [1..l]]; // total + m...(m-l+1) --- & fails if range empty! end while; end for; print m, n, special; // should sum to A180238 /* Output Time: 0.000 1 15 [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ] Time: 0.380 2 15 [ 1, 2, 4, 8, 14, 24, 36, 54, 76, 104, 136, 178, 224, 282, 346, 418 ] Time: 7.770 3 15 [ 1, 3, 9, 27, 75, 189, 447, 951, 1911, 3621, 6513, 11103, 18267, 29013, 44691, 67251 ] 3 15 [ [ 0, 0, 0, 6, 24, 78, 186, 372, 876, 1632, 3024, 5310, 8496, 13344, 21186, 0 ], [ 0, 0, 6, 18, 36, 78, 150, 306, 420, 792, 1338, 2082, 3228, 4830, 7050, 0 ], [ 1, 3, 3, 3, 9, 15, 33, 63, 153, 219, 261, 351, 585, 879, 933, 0 ] ] */