print(` Version of May 20 1999`):
print(` Previous Version was written Dec. 4, 1997`):
print(`This is DAVID_IAN, A Maple package`):
print(`that implements the Goulden-Jackson Cluster method.`):
print(`Written by John Noonan and Doron Zeilberger.`):
lprint(``):
print(`It is one of the packages accompanying Noonan and Zeilberger's `):
print(`article: The Goulden-Jackson Cluster Method: Extensions,`):
print(`Applications, and Implementations, where the method`):
print(`is explained in detail, and then extended and generalized in`):
print(`various directions. `):
lprint(``):
print(`The present package only deals with the original method.`):
print(`Other packages deal with the extensions`):
lprint(``):
print(`Ian Goulden and David Jackson's original treatment can be found `):
print(`in their paper that appeared in J. London Math. Soc.(2)20(1979),`):
print(` 567-576. It is also described in their book`):
print(`Combinatorial Enumeration, John  Wiley, 1983, New York`):
print(`See also R. Stanley's Enumerative Combinatorics, v.1, 266, 280.`):
print(`and Guibas, Odlyzko,JCT(A)30(1981),183-208.`):
lprint(``):
print(`The paper, and the packages, including the present one`):
print(`is available from http://www.math.temple.edu/~zeilberg`):
lprint(``):
print(`For general help, and a list of the available functions,`):
print(` type ezra();. For specific help type ezra(procedure_name) `):
lprint(``):
 
 
ezra:=proc()
if args=NULL then
print(`DAVID_IAN`):
print(`A Maple package that implements the Goulden-Jacskon Cluster method`):
print(` for finding generating functions and series expansions`):
print(`for the number of words avoiding a prescribed, finite, set of`):
print(`mistakes. Contains Procedures , BestLastPlay `):
print(` EmpiricalGJseries `):
print(` GJseries,GJgf,GJs, GJst, GJgfDetail, GJstDetail, Penney. `):
print(` PenneyGames, SetOfCorrectWords, wGJseries `):
print(` For specific help type ezra(procedure_name) `):
fi:
 
if nops([args])=1 and op(1,[args])=EmpiricalGJseries then
print(`EmpiricalGJseries(alphabet,MISTAKES1,L): Does what`):
print(`GJseries(nops(alphabet),MISTAKES,L) does, but`):
print(`empirically. Don't make L to large, unless there `):
print(`are lots of mistakes. It is meant to be a check`):
print(`on the validity of the (usually) much much faster`):
print(`and should not be used except for curiosity`):
fi:
 
 
if nops([args])=1 and op(1,[args])=SetOfCorrectWords then
print(`SetOfCorrectWords(Alphabet,MISTAKES,n): given an alphabet, Alphabet`):
print(`and a set of mistakes, MISTAKES, finds, empirically, the`):
print(`set of words that that do not have factors that belong to MISTAKES`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`wGJseries` then
print(`wGJseries(alphabet,MISTAKES1,x,L): The weighted series expansion`):
print(`of the words of length n in the alphabet alphabet, avoiding`):
print(`the set of mistakes MISTAKES1, up to L terms, where the`):
print(`i_th term is the weight-enumerator`):
fi:
 
if nops([args])=1 and op(1,[args])=`BestLastPlay` then
print(`BestLastPlay(alphabet,FirstMoves,Probs,R)`):
print(`In Penney Ante, given a die whose sides are labelled by letters`):
print(`of the alphabet alphabet, loaded like in the list`):
print(`Probs, and a list of moves of the previous players`):
print(`FirstMoves, returns the best`):
print(`move with R letters`):
print(`followed by the list of probabilities of winning`):
print(`of all the players`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`PenneyGames` then
print(`PenneyGames(alphabet,MISTAKES,Probs,K): playes K Penney-Ante`):
print(`games, and outputs the scores of each respective winner`):
print(`according to the order of MISTAKES`):
fi:
 
if nops([args])=1 and op(1,[args])=`Penney` then
print(`Penny(alphabet,WORDS,Probs): computes the relative `):
print(`probabilites of winning a Penny-Ante game where each`):
print(`player picks a word beforehand, and a letter is chosen`):
print(`at random until one of the words shows up`):
print(`alphabet is the alphabet (usually [H,T]), given in terms of`):
print(` a list, WORDS are the words, and Probs is the list of`):
print(`probabilites that each of the letters of the alphabet`):
print(`has of showing up.`):
print(`The output gives the respecive probablities of the words`):
print(`of WORDS of showing up first`):
lprint(``):
print(`For example, to do the Penney-Ante game described in`):
print(` Concrete Math do Penney([H,T],[[H,H,T],[H,T,T]],[1/2,1/2]):`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`GJseries` then
print(`GJseries(NUMLETTERS,MISTAKES,L) finds the first L terms `):
print(`of the sequence: Number of words of length n, having no factors`):
print(`(i.e. subsequences of consecutive letters) that belong to the`):
print(`set of lists MISTAKES, where the number of letters of the alphabet`):
print(` is NUMLETTERS. For example to find the sequence a_n (0<=n<=15)`):
print(`where a_n:= number of words in {1,2,3} avoiding 123,231,312`):
print(` type GJseries(3,{[1,2,3],[2,3,1],[3,1,2]},15)`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`GJgf` then
print(`GJgf(alphabet,MISTAKES,x,t) finds the generating function`):
print(`F(x[l_1],x[l_2],...,x[l_r];t), the weight enumerator of `):
print(`of all words in the alphabet=[l_1,..,l_r], where the`):
print(`the weight of a word is the product of x[its letters] times`):
print(`t^(#number of factors that belong to the list MISTAKES`):
print(`the set of mistakes, MISTAKES should be such that if it contains`):
print(` a word, it cannot contain any of its proper subwords`):
print(`For example to get the g.f. for the set of words`):
print(`in the alphabet {1,2,3}, avoiding 123,231,312, type`):
print(` type GJgf({1,2,3},{[1,2,3],[2,3,1],[3,1,2]},x,t)`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`GJgfDetail` then
print(`GJgfDetail(alphabet,MISTAKES,x,t) finds the generating function`):
print(`F(x[l_1],x[l_2],...,x[l_r];t), the weight enumerator of `):
print(`of all words in the alphabet=[l_1,..,l_r], where the`):
print(`the weight of a word is the product of x[its letters] times`):
print(`product of t[w] for all factors w that belong to MISTAKES.`):
print(`the set of mistakes, MISTAKES should be such that if it contains`):
print(` a word, it cannot contain any of its proper subwords`):
print(`For example to get the g.f. for the set of words`):
print(`in the alphabet {1,2,3}, avoiding 123,231,312, type`):
print(` type GJgfDetail({1,2,3},{[1,2,3],[2,3,1],[3,1,2]},x,t)`):
fi:
 
 
 
if nops([args])=1 and op(1,[args])=`GJst` then
print(`GJst(alphabet,MISTAKES,s,t) finds the generating function`):
print(`F(s,t), the weight enumerator of `):
print(`of all words in the alphabet=[l_1,..,l_r], where the`):
print(`the weight of a word is s^(length) times `):
print(`t^(#number of factors that belong to the list MISTAKES`):
print(`the set of mistakes, MISTAKES should be such that if it contains`):
print(` a word, it cannot contain any of its proper subwords`):
print(`For example to get the g.f. for the set of words`):
print(`in the alphabet {1,2,3}, avoiding 123,231,312, type`):
print(` type GJst({1,2,3} ,{[1,2,3],[2,3,1],[3,1,2]},s,t)`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`GJstDetail` then
print(`GJstDetail(alphabet,MISTAKES,s,t) finds the generating function`):
print(`F(s,t[mistakes]), the weight enumerator of `):
print(`of all words in the alphabet=[l_1,..,l_r], where the`):
print(`the weight of a word is s^(length) times `):
print(`product of t[w] for all factors w that belong to MISTAKES.`):
print(`the set of mistakes, MISTAKES should be such that if it contains`):
print(` a word, it cannot contain any of its proper factors.`):
lprint(``):
print(`For example to get the g.f. for the set of words`):
print(`in the alphabet {1,2,3}, avoiding 123,231,312, type`):
print(` type GJstDetail({1,2,3} ,{[1,2,3],[2,3,1],[3,1,2]},s,t)`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`GJs` then
print(`GJs(alphabet,MISTAKES,s) finds the generating function`):
print(`g(s), is the g.f. sum(a[i]*s^i,i=0..infty)  `):
print(`where a[i]:=number of words of length i without any mistakes`):
print(`from MISTAKES`):
print(`For example to get the g.f. for a_n:=number of words of length n`):
print(`in the alphabet {1,2,3}, avoiding 123,231,312, type`):
print(` type GJs({1,2,3},{[1,2,3],[2,3,1],[3,1,2]},s)`):
fi:
 
 
 
 
 
 
end:
 
 
 
 
#A program to find the gen. fun. of all words in an alphabet
#given in a list alphabet, according to the number of mistakes
#where the list of mistakes is MISTAKES, a list of lists
#where each element is given in as a list. It is based on
#Goulden-Jackson-Zeilberger, see Stanley's book, 
#Enumerative combinatorics I, ex. 14.a in
#chapter 4, pp. 266-267
 
 
#overlap is a procedure that given two words u and v
#computes the weight-enumerator of all v\suffix(u),
#for all suffixes of u that are prefixes of v.
overlap:=proc(u,v,x)
local i,j,k,lu,gug:
 
lu:=0:
 
for i from 2 to nops(u) do
  for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v))
         do
   :
  od:   
 
if j-i=nops(v) and u<>v then
 ERROR(v,`is a subword of`,u,`illegal input`):
fi:
 
 
 
 if j=nops(u)+1 and (i>1 or j>2) then
  gug:=1:
 
  for k from j-i+1 to nops(v) do
   gug:=gug*x[op(k,v)]:
  od:
 
  lu:=lu+gug:
 
fi:
 
od:
 
lu:
 
end:
 
#overlapz is a procedure that given two words u and v, and a variable s
#computes the weight-enumerator of all v\suffix(u),
#for all suffixes of u that are prefixes of v, but with uniform weight s
overlapz:=proc(u,v,s)
local i,j,lu,gug:
 
lu:=0:
 
for i from 2 to nops(u) do
  for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v))
         do
   :
  od:   
 
if j-i=nops(v) and u<>v then
 ERROR(v,`is a subword of`,u,`illegal input`):
 fi:
 
 
 
 if j=nops(u)+1 and (i>1 or j>2) then
  gug:=1:
 
   gug:=gug*s^(nops(v)-(j-i)):
 
  lu:=lu+gug:
 
fi:
 
od:
 
lu:
 
end:
 
 
#findeq sets up the equ C[v]= x[v]+t*Sum_u overlap(u,v,x) *C[u]
 
findeq:=proc(v,MISTAKES,C,t,x)
local eq,i,u:
eq:=t:
 
for i from 1 to nops(v) do
 eq:=eq*x[op(i,v)]:
od:
 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq+t*overlap(u,v,x)*C[op(u)]:
od:
 
C[op(v)]-eq=0:
 
end:
 
 
#findeqz sets up the equ C[v]= s+t*Sum_u overlap(u,v,x) *C[u]
 
findeqz:=proc(v,MISTAKES,C,s)
local eq,i,u:
eq:=-1:
 
for i from 1 to nops(v) do
 eq:=eq*s:
od:
 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq-overlapz(u,v,s)*C[op(u)]:
od:
 
C[op(v)]-eq=0:
 
end:
 
 
GJgf:=proc(alphabet,MISTAKES,x,t)
local v,eq, var,i,lu,C:
 
if Hakten(MISTAKES)<>MISTAKES then
ERROR(`The set of mistakes is not minimal, use another package`):
fi:
 
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeq(v,MISTAKES,C,t,x)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
lu:=1:
 
for i from 1 to nops(alphabet) do
 lu:=lu-x[op(i,alphabet)]:
od:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
od:
 
lu:=normal(1/lu):
subs(t=t-1,lu): 
end:
 
  
GJst:=proc(alphabet,MISTAKES,s,t)
local v,eq, var,i,lu,C,x:
 
if Hakten(MISTAKES)<>MISTAKES then
ERROR(`The set of mistakes is not minimal, use another package`):
fi:
 
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeq(v,MISTAKES,C,t,x)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
lu:=1:
 
for i from 1 to nops(alphabet) do
 lu:=lu-s:
od:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
od:
 
for i from 1 to nops(alphabet) do
lu:=subs(x[op(i,alphabet)]=s,lu):
od: 
lu:=normal(1/lu):
normal(subs(t=t-1,lu)): 
end:
 
 
 
 
GJs:=proc(alphabet,MISTAKES1,s)
local v,eq, var,i,lu,C,MISTAKES:
 
MISTAKES:=Hakten(MISTAKES1):
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeqz(v,MISTAKES,C,s)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
lu:=1:
 
for i from 1 to nops(alphabet) do
 lu:=lu-s:
od:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
od:
 
normal(1/lu):
 
end:
 
 
 
 
 
#given a word, computes all its shifts
 
shift1:=proc(u)
local gu,i:
gu:=[u]:
 
for i from 2 to nops(u) do
 gu:=[op(gu),[op(i..nops(u),u),op(1..i-1,u)]]:
od:
 
gu:
 
end:
 
 
#overlap1 finds the list of places i in v s.t. v[1],...,v[i] is a suffix of u
overlap1:=proc(u,v)
local i,gu:
gu:=[]:
 
for i from 1 to nops(v)-1 do
 if nops(u)-i+1>1 and [op(1..i,v)]=[op(nops(u)-i+1..nops(u),u)] then
  gu:=[op(gu),nops(v)-i]:
 fi:
od:
 
gu:
 
end:
 
 
 
#Given a set of mistakes overlapmat1 is a list of lists A such that 
#op(i,op(j,A)) is the overlap(op(j,MISTAKES),op(i,MISTAKES)))
 
overlapmat1:=proc(MISTAKES)
local gu,u,v,lu,i,j:
gu:=[]:
 
for i from 1 to nops(MISTAKES) do
v:=op(i,MISTAKES):
lu:=[]:
 
  for j from 1 to nops(MISTAKES) do
   u:=op(j,MISTAKES):
    lu:=[op(lu),overlap1(u,v)]:
  od:
 
 gu:=[op(gu),lu]:
od:
 
gu:
 
end:
 
 
 
GJseries:=proc(NUMLETTERS,MISTAKES1,L)
local A, MAXOVER,i,j,k,CU,NUM,kuku,kuku1,lu,pip,CUU,KHADAS,gagi,r,v,TOT,TOV,
i1,khad,MISTAKES:
 
MISTAKES:=Hakten(MISTAKES1):
if min(seq(nops(MISTAKES[i1]),i1=1..nops(MISTAKES)))<=1 then
ERROR(`The Mistakes should all be of length at least 2`):
fi:

TOT:=[]:
 
NUM:=nops(MISTAKES):
 
A:=overlapmat1(MISTAKES);
 
MAXOVER:=0:
 
 
for i from 1 to NUM do
  kuku:=op(i,A):
 for j from 1 to NUM do
  kuku1:=op(j,kuku):
       if nops(kuku1)>0 and max(op(kuku1))>MAXOVER then
      MAXOVER:=max(op(kuku1)):
    fi:
 od:
od:
 
#initialize CU
 
 
 CU:=[]:
for i from 1 to MAXOVER do
 lu:=[]:
 
 for j from 1 to NUM do
  lu:=[op(lu),0]:
 od:
   CU:=[op(CU),lu]:
od:
 
 
for k from 1 to L do
 KHADAS:=[]:
 
  for i from 1 to NUM do
   v:=op(i,MISTAKES):
 
    if nops(v)=k then
     pip:= -1:
    else
       pip:=0:
    fi:
 
 for j from 1 to NUM do
   gagi:=op(j,op(i,A)):
  for r from 1 to nops(gagi) do
    pip:=pip-op(j,op(op(r,gagi),CU))
  od:
od:
 
KHADAS:=[op(KHADAS),pip]:
 
od:
 
TOT:=[op(TOT),convert(KHADAS,`+`)]:
CUU:=[KHADAS]:
 
for i from 1 to MAXOVER-1 do
 CUU:=[op(CUU),op(i,CU)]:
od:
 
CU:=CUU:
od:
 
 
TOV:=[1]:
 
for r from 1 to L do
 
khad:=NUMLETTERS*op(nops(TOV),TOV):
 for i1 from 2 to r do
  khad:=khad+op(r-i1+1,TOV)*op(i1,TOT):
 od:
 
TOV:=[op(TOV),khad]:
 
od:
 
TOV:
 
end:
 
 
 
 
findeqDetail:=proc(v,MISTAKES,C,t,x)
local eq,i,u:
eq:=t[op(v)]:
 
for i from 1 to nops(v) do
 eq:=eq*x[op(i,v)]:
od:
 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq+t[op(v)]*overlap(u,v,x)*C[op(u)]:
od:
 
C[op(v)]-eq=0:
 
end:
 
 
GJgfDetail:=proc(alphabet,MISTAKES,x,t)
local v,eq, var,i,lu,C:
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeqDetail(v,MISTAKES,C,t,x)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
lu:=1:
 
for i from 1 to nops(alphabet) do
 lu:=lu-x[op(i,alphabet)]:
od:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
od:
 
lu:=normal(1/lu):
for i from 1 to nops(MISTAKES) do
v:=op(i,MISTAKES):
lu:=subs(t[op(v)]=t[op(v)]-1,lu): 
od:
normal(lu):
end:
 
 
findeqzDetail:=proc(v,MISTAKES,C,t,s)
local eq,i,u:
eq:=t[op(v)]*s^(nops(v)):
 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq+t[op(v)]*overlapz(u,v,s)*C[op(u)]:
od:
 
C[op(v)]-eq=0:
 
end:
 
 
GJstDetail:=proc(alphabet,MISTAKES,s,t)
local v,eq, var,i,lu,C:
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeqzDetail(v,MISTAKES,C,t,s)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
lu:=1-s*nops(alphabet):
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
od:
 
lu:=normal(1/lu):
for i from 1 to nops(MISTAKES) do
v:=op(i,MISTAKES):
lu:=subs(t[op(v)]=t[op(v)]-1,lu): 
od:
normal(lu):
end:
 
 
Penney:=proc(alphabet,MISTAKES,Probs)
local v,eq, var,i,C,x,t,ANSWER,mu,ANSWER1:
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {subs(t=-1,findeq(v,MISTAKES,C,t,x))}:
 var:=var union {C[op(v)]}:
od:
 
for i from 1 to nops(alphabet) do
eq:=subs(x[op(i,alphabet)]=op(i,Probs),eq):
od:
 
 
var:=solve(eq,var):
 
if var=NULL then
ERROR(`Something is fishy`):
fi:
 
 
ANSWER:=[]:
 
for i from 1 to nops(MISTAKES) do
v:=op(i,MISTAKES):
ANSWER:=[op(ANSWER),subs(var,C[op(v)])]:
od:
 
mu:=0:
 
for i from 1 to nops(ANSWER) do
mu:=mu+op(i,ANSWER):
od:
 
ANSWER1:=[]:
 
 
for i from 1 to nops(ANSWER) do
ANSWER1:=[op(ANSWER1),normal(op(i,ANSWER)/mu)]:
od:
ANSWER1:
end:
 
 
 
#issubword(u,v) returns 1 if v is a subword of u, otherwise 0
 
issubword:=proc(u,v)
local i,j:
 
 
for i from 1 to nops(u) do
  for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v))
         do
   :
  od:   
 
if j-i=nops(v)  then
 RETURN(1):
 fi:
od:
 
0:
 
end:
 
#superflous(B,w) given a set of bad words, and one of its
#members, finds whether it is superflous,
#that is, whether it is a factor of another one, and deletes the larger one
 
superflous:=proc(B,w)
local i,v:
 
if not type(B,set) or not type(w,list) then
 ERROR(`Bad input`):
fi:
 
if not member(w,B) then
 ERROR(`Second argument must belong to first arg.`):
fi:
 
for i from 1 to nops(B) do
v:=op(i,B):
if v<>w and issubword(w,v)=1 then
 RETURN(1):
fi:
od:
0:
end:
 
hakten:=proc(B)
local w,i:
 
for i from 1 to nops(B) do
w:=op(i,B):
 
if superflous(B,w)=1 then
 RETURN(B minus {w}):
fi:
od:
B:
end:
 
#Hakten(B) removes all superflous words
 
Hakten:=proc(B)
local B1,B2:
 
B1:=B:
 
B2:=hakten(B1):
 
while B1<>B2 do
B1:=B2:
B2:=hakten(B2):
od:
 
B2:
end:
 
 
#PenneyGames(alphabet,MISTAKES,Probs,K): playes K Penney-Ante
#games, and outputs the scores of each respective winner
#according to the order of MISTAKES:
 
PenneyGames:=proc(alphabet,MISTAKES,Probs,K)
local lu,i,winner1:
 
lu:=[]:
for i from 1 to nops(MISTAKES) do
 lu:=[op(lu),0]:
od:
 
for i from 1 to K do
winner1:=PenneyGame1(alphabet,MISTAKES,Probs):
lu:=[op(1..winner1-1,lu),op(winner1,lu)+1,op(winner1+1..nops(lu),lu)]:
od:
lu:
end:
 
 
#PenneyGame1(alphabet,MISTAKES,Probs,K): performs 1 Penne-Ante
#game where the die is labelled by the letters of the set
#alphabet, the probablities are given by Probs
 
PenneyGame1:=proc(alphabet,MISTAKES,Probs)
local w,ra,j,lu,thr:
ra:=ldie(Probs):
lu:=ra[2]:
ra:=ra[1]:
 
w:=[]:
 
while nigmar(w,MISTAKES)=0 do
thr:=ra():
 
for j from 1 to nops(lu) do
if thr<=op(j,lu) then
 w:=[op(w),op(j,alphabet)]:
 break:
fi:
od:
 
od:
nigmar(w,MISTAKES):
end:
 
#loaded die given by Probs
 
 
ldie:=proc(Probs)
local ku,Probs1,i,mu,lu,ra:
 
if nops(Probs)<2 then
ERROR(`Input is a list with at least two elements`):
fi:
 
 
mu:=0:
 
for i from 1 to nops(Probs) do
if op(i,Probs)<0 then
ERROR(`Probabilities must be non-negative`):
fi:
mu:=mu+op(i,Probs):
od:
 
Probs1:=[]:
 
for i from 1 to nops(Probs) do
Probs1:=[op(Probs1),op(i,Probs)/mu]:
od:
 
mu:=lcm(denom(op(1,Probs1)),denom(op(2,Probs1))):
for i from 3 to nops(Probs1) do
mu:=lcm(mu,denom(op(i,Probs1))):
od:
 
 
 
lu:=numer(mu*op(1,Probs1)):
 
ku:=[]:
for i from 1 to nops(Probs1)-1 do
ku:=[op(ku),lu]:
lu:=lu+numer(mu*op(i+1,Probs1)):
od:
ku:=[op(ku),lu]:
 
ra:=rand(1..mu):
ra,ku:
 
end:
 
 
 
 
#for i from 1 to nops(lu) do
#if ra<=op(i,lu) then
#RETURN(i):
#fi:
#od:
 
 
#nigmar(w,MISTAKES): given a word w and a set of words MISTAKES
#decides whether the end of w coincides with one of the words
#of MISTAKES and then outputs which one, otherwise gives 0
 
nigmar:=proc(w,MISTAKES)
local i,mila:
 
for i from 1 to nops(MISTAKES) do
mila:=op(i,MISTAKES):
 
if nops(w)>=nops(mila) and [op(nops(w)-nops(mila)+1..nops(w),w)]=mila then
RETURN(i):
fi:
od:
0:
end:
 
 
#Seqk(alphabet,k): Given a set of letters and a non-negative
#integer k, outputs the set of words of length k in that alphabet
Seqk:=proc(alphabet,k)
local lu,j,i,w,gu:
 
if k=0 then
 RETURN({[]}):
fi:
lu:=Seqk(alphabet,k-1):
 
gu:={}:
 
for i from 1 to nops(lu) do
w:=op(i,lu):
for j from 1 to nops(alphabet) do
gu:=gu union {[op(w),op(j,alphabet)]}:
od:
od:
gu:
end:
 
 
 
#BestLastPlay(alphabet,FirstMoves,Probs,R)
#In Penney Ante, given a die whose sides are labelled by letters
#of the alphabet alphabet, loaded like in the list
#Probs, and a list of moves of the previous players
#FirstMoves, returns the best
#move with R letters, followed by the list of probabilities of winning
#of all the players
 
 
BestLastPlay:=proc(alphabet,FirstMoves,Probs,R)
local MISTAKES,i,lu,w,ma,gu,aluf,histab:
 
ma:=0:
 
lu:=Seqk(alphabet,R) minus convert(FirstMoves,set):
 
for i from 1 to nops(lu) do
w:=op(i,lu):
MISTAKES:=[op(FirstMoves),w]:
 
 
if Hakten(convert(MISTAKES,set))=convert(MISTAKES,set) then
gu:=Penney(alphabet,MISTAKES,Probs):
 
if op(nops(MISTAKES),gu)>ma then
 ma:=op(nops(MISTAKES),gu):
 aluf:=w:
 histab:=gu:
fi:
fi:
od:
 
aluf,histab:
 
end:
 
#######################NEW STUFF
 
 
mishkal:=proc(w,x)
local i,gu:
gu:=1:
 
for i from 1 to nops(w) do
gu:=gu*x[w[i]]:
od:
gu:
end:
 
 
#woverlap1 finds the list of places i in v s.t. 
#v[1],...,v[i] is a suffix of u accompanied by the 
#weight of the leftover
 
woverlap1:=proc(u,v,x)
local i,gu:
gu:=[]:
 
for i from 1 to nops(v)-1 do
 if nops(u)-i+1>1 and [op(1..i,v)]=[op(nops(u)-i+1..nops(u),u)] then
  gu:=[op(gu),[nops(v)-i,mishkal([op(i+1..nops(v),v)],x)]]:
 fi:
od:
 
gu:
 
end:
 
 
#Given a set of mistakes, woverlapmat1(MISTAKES,x)
# is a list of lists A such that 
#op(i,op(j,A)) is the (the weighted overlap)
#woverlap(op(j,MISTAKES),op(i,MISTAKES)))
 
woverlapmat1:=proc(MISTAKES,x)
local gu,u,v,lu,i,j:
gu:=[]:
 
for i from 1 to nops(MISTAKES) do
v:=op(i,MISTAKES):
lu:=[]:
 
  for j from 1 to nops(MISTAKES) do
   u:=op(j,MISTAKES):
    lu:=[op(lu),woverlap1(u,v,x)]:
  od:
 
 gu:=[op(gu),lu]:
od:
 
gu:
 
end:
 
 
 
#wGJseries(alphabet,MISTAKES1,x,L): The weighted series expansion
#of the words of length n in the alphabet alphabet, avoiding
#the set of mistakes MISTAKES1, up to L terms, where the
#i_th term is the weight-enumerator
 
wGJseries:=proc(alphabet,MISTAKES1,x,L)
local A, MAXOVER,i,j,k,CU,NUM,kuku,kuku1,lu,pip,CUU,KHADAS,gagi,r,v,TOT,TOV,
i1,khad,MISTAKES,lok,SUMLETTERS,kuku1n:
 
SUMLETTERS:=0:
for i from 1 to nops(alphabet) do
SUMLETTERS:=SUMLETTERS+x[op(i,alphabet)]:
od:
 
MISTAKES:=Hakten(MISTAKES1):
TOT:=[]:
 
NUM:=nops(MISTAKES):
 
A:=woverlapmat1(MISTAKES,x);
 
MAXOVER:=0:
 
 
for i from 1 to NUM do
  kuku:=op(i,A):
 for j from 1 to NUM do
  kuku1:=op(j,kuku):
 
 kuku1n:={}:
 
 for lok from 1 to nops(kuku1) do
  kuku1n:=kuku1n union {op(lok,kuku1)[1]}:
 od:
 
       if nops(kuku1)>0 and max(op(kuku1n))>MAXOVER then
      MAXOVER:=max(op(kuku1n)):
    fi:
 od:
od:
 
#initialize CU
 
 
 CU:=[]:
for i from 1 to MAXOVER do
 lu:=[]:
 
 for j from 1 to NUM do
  lu:=[op(lu),0]:
 od:
   CU:=[op(CU),lu]:
od:
 
 
for k from 1 to L do
 KHADAS:=[]:
 
  for i from 1 to NUM do
   v:=op(i,MISTAKES):
 
    if nops(v)=k then
     pip:= -1:
      for lok from 1 to nops(v) do
       pip:=pip*x[v[lok]]:
      od:
 
    else
       pip:=0:
    fi:
 
 for j from 1 to NUM do
   gagi:=op(j,op(i,A)):
  for r from 1 to nops(gagi) do
    pip:=pip-gagi[r][2]*op(j,op(op(r,gagi)[1],CU))
  od:
od:
 
KHADAS:=[op(KHADAS),expand(pip)]:
 
od:
 
TOT:=[op(TOT),convert(KHADAS,`+`)]:
CUU:=[KHADAS]:
 
for i from 1 to MAXOVER-1 do
 CUU:=[op(CUU),op(i,CU)]:
od:
 
CU:=CUU:
od:
 
 
TOV:=[1]:
 
for r from 1 to L do
 
khad:=SUMLETTERS*op(nops(TOV),TOV):
 for i1 from 2 to r do
  khad:=khad+op(r-i1+1,TOV)*op(i1,TOT):
 od:
 
TOV:=[op(TOV),expand(khad)]:
 
od:
 
TOV:
 
end:
 
 
#####
 
#isbad(mila,MISTAKES): returns true if one of the suffices of
#mila belongs to the set of MISTAKES
isbad:=proc(mila,MISTAKES)
local i,shegi:
 
for i from 1 to nops(MISTAKES) do
shegi:=op(i,MISTAKES):
 
if nops(mila)>=nops(shegi) and 
 [op(nops(mila)-nops(shegi)+1..nops(mila),mila)]=shegi then
 RETURN(true):
fi:
 
od:
 
false:
 
end:
 
 
 
#SetOfCorrectWords(Alphabet,MISTAKES,n): given an alphabet, Alphabet
#and a set of mistakes, MISTAKES, finds, empirically, the
#set of words that that do not have factors that belong to MISTAKES
SetOfCorrectWords:=proc(Alphabet,MISTAKES,n)
local MISTAKES1,gu,mu,i,j,cand:
option remember:
MISTAKES1:=Hakten(MISTAKES):
 
if n=0 then
RETURN({[]}):
fi:
 
mu:=SetOfCorrectWords(Alphabet,MISTAKES,n-1):
 
gu:={}:
 
for i from 1 to nops(mu) do
 for j from 1 to nops(Alphabet) do
   cand:=[op(op(i,mu)),op(j,Alphabet)]:
 
     if not isbad(cand,MISTAKES1) then
          gu:=gu union {cand}:
     fi:
 od:
od:
 
gu:
 
end:
 
 
 
#EmpiricalGJSeries(alphabet,MISTAKES1,L): Does what
#GJseries(nops(alphabet),MISTAKES,L) does, but
#empirically. Don't make L to large, unless there 
#are lots of mistakes. It is meant to be a check
#on the validity of the (usually) much much faster
#and should not be used except for curiosity
EmpiricalGJseries:=proc(alphabet,MISTAKES,L)
local n,gu:
gu:=[]:
 
for n from 0 to L do
gu:=[op(gu),nops(SetOfCorrectWords(alphabet,MISTAKES,n))]:
od:
gu:
end: