######################################################################
##Litt.txt: Save this file as  Litt.txt                              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read Litt.txt<Enter>                               #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: May 2024
 
print(`Created:  May 2024`):
print(`Previous version: June 10, 2024, adding procedure Mihai, implementing  the amazing  formula of Mihai Nica`):
print(`This version: June 14, 2024, adding procedure HO(m,word) and NicaConj(A,B)`):
print(` This is Litt.txt `):
print(`It is a Maple package that accompanies the article `):
print(`How to Answer Questions of the Type: If you toss a coin n times, how likely is HH to show up more than HT?`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):

print(`available from http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/litt.html`):

print(``):
print(`Please report bugs to DoronZeil@gmail.com`):
print(``):
 print(`The most current version of this  package `):
 print(` is  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/Litt.txt  .`):

print(`-------------------------------`):
 print(`For a list of the main procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-------------------------------`):


 print(`For a list of the  other procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(`-------------------------------`):

 print(`For a list of the Checking procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

print(`-------------------------------`):

with(combinat):

ezraC:=proc()

if args=NULL then
 print(` The  checking procedures are:   ABseqBF, CheckGFwtE,  TestMihali, TestMihai1, TestNicaConj1, TestNicaConj`):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AvsB,  BETc,  Book1, CoeffXn, DOaz, EWords, findrec, Findrec, GFwtE,  GFwtEold, HO, MomsAB, NuSS, NuSS1, SeqFromRec, TestConjVar, WhoWonVar, Words, WtE  `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: ABseq, Book, CounterBets, Info, LD, Mihai, NicaBook, RECaz, WhoWon,  WhoWonN, WhoWonV, WW `):
 print(` `):



elif nops([args])=1 and op(1,[args])=ABseq then
print(`ABseq(n,m,A,B): the first n terms of the sequence (starting at n1=1) of the sequence AvsB(n1,m,A,B), the clever way. Try:`):
print(`ABseq(10,2,{[1,1]},{[1,2]});`):

elif nops([args])=1 and op(1,[args])=ABseqBF then
print(`ABseqBF(n,m,A,B): the first n terms of the sequence (starting at n1=1) of the sequence AvsB(n1,m,A,B) by brute force. For checking only. Try:`):
print(`ABseqBF(10,2,{[1,1]},{[1,2]});`):

elif nops([args])=1 and op(1,[args])=AvsB then
print(`AvsB(n,m,A,B): inputs pos. integers n and m, and set of words A and B in the alphabet {1,...,m} finds the number of n-letter words in {1,..,m} where`):
print(`there are more occurrences of members of A then B. Try:`):
print(`AvsB(5,2,{[1,1]},{[1,2]});`):


elif nops([args])=1 and op(1,[args])=BETc then
print(`BETc(BET,m): All the equivalence classes of B of words of length nops(BET) for B under "having the same GFwtE(m,{BET},{B},a,b,z). Try`):
print(`BETc([1,1,1],2);`):

elif nops([args])=1 and op(1,[args])=Book then
print(`Book(m,k,K): A book about the relative probablity of winning Daniel Litt bets for all k-letter words in the alphabet {1,...,m} using K terms. Try:`):
print(`Book(2,3,5000);`):

elif nops([args])=1 and op(1,[args])=Book1 then
print(`Book1(m,BET,K): inputs a positive integer m (at least 2), and, BET, ( a word in the alphabet {1,...,m} of length at least 2)`):
print(`and a large positive intetger K (at least 1000) outputs`):
print(`a book of the winning chances in the Litt game of all k-letter words in {1,...,m} versus BET. Try:`):
print(`Book1(2,[1,1,1,1],5000);`):

elif nops([args])=1 and op(1,[args])=CheckGFwtE then
print(`CheckGFwtE(m,A,B,n): checks GFwtE(m,A,B,a,b,x) up to the n-th coefficient. Try:`):
print(`CheckGFwtE(2,{[1,1]},{[1,2]},10);`):

elif nops([args])=1 and op(1,[args])=CoeffXn then
print(`CoeffXn(f,x,n): given a rational function whose denominator is a power of (1-x) finds the coeff. of x^n in its Maclaurin expansion. Try:`):
print(`CoeffXn(x^4/(1-x)^3,x,n);`):

elif nops([args])=1 and op(1,[args])=CounterBets then
print(`CounterBets(m,A): inputs a word A in {1,..,m} and outputs the randked list  of the counterbets, from best to worst. Try:`):
print(`CounterBets(2,[1,1,1]);`):

elif nops([args])=1 and op(1,[args])=DOaz then
print(`DOaz(m,A,B,z,Dz): Given a positive integer`):
print(`m and two sets of words A and B of words in {1,...,m} of the same length in the`):
print(`uses the continuous Almkvist-Zeilberger algorithm,`):
print(`to find a linear differential operator with`):
print(`polynomial coefficients, in terms of z and the`):
print(`differentiation w.r.t. z, Dz, that annihilates`):
print(`the (ordinary) generating function`):
print(`of the sequence`):
print(`"number of words in {1, ...m} of length n that have more occurrences of `):
print(`words in A then words in B`):
print(`For example for the original Daniel Litt problem`):
print(`DOaz(2,{[1,1]},{[1,2]},t,Dt);`):


elif nops([args])=1 and op(1,[args])=EWords then
print(`EWords(n,m): A set of representatives of  n-letter words in {1,...,m}. Try: `):
print(`EWords(4,3);`):

elif nops([args])=1 and op(1,[args])=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nops([args])=1 and op(1,[args])=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs.`):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nops([args])=1 and op(1,[args])=GFwtE then
print(`GFwtE(m,A,B,a,b,x): inputs a positive integer m, sets of words in {1,...,m}, A,B, variables, a and b, and another variable x, finds the rational function`):
print(`Sum(WtE(n,m,A,B,a,b)*x^n,n=0..infinity); It use the Goulden-Jackson cluster method. Try:`):
print(`GFwtE(2,{[1,1]},{[1,2]},a,b,x);`):


elif nops([args])=1 and op(1,[args])=GFwtEold then
print(`GFwtEold(m,A,B,a,b,x): Same as GFwtE(m,A,B,a,n,x), but slower, using the direct. approach. Try:`):
print(`GFwtEold(2,{[1,1]},{[1,2]},a,b,x);`):


elif nops([args])=1 and op(1,[args])=HO then
print(`HO(m,w): the Hanging Overlap of the word w according to Mihai Nica. Try:`):
print(`HO(2,[1,2,1,2,1]);`):

elif nops([args])=1 and op(1,[args])=Info then
print(`Info(m,A,B,r,K,n,N,C): Given a positive integer m>=2, and two sets of words (all consisting of the same length) in the alphabet {1,...,m}, A,B, and a large positive integer K and a positive integer`):
print(`r, outputs the pair of lists`):
print(`[[ope1,Ini1,Asy1],[ope2,Ini2,Asy2]]. `):
print(`Where ope1 is the operator (of complexity at most C) annihilating the sequence 2^(n-1)-Number of n-letter words with more A's than B, Ini1 the initial condition and Asy1 the`):
print(`asympotics to order r, `):
print(`[ope2,Ini2,Asy2] same but with 2^(n-1)-Number of n-letter words with more B's than A. Try:`):
print(`Info(2,{[1,1]},{[1,2]},5,5000,n,N,10);`):



elif nops([args])=1 and op(1,[args])=LD then
print(`LD(A,B,N): inputs two sets of words A and B consisting of words all of the same`):
print(`length, and a positive integer N outputs the list of `):
print(`length N such that if Pr(H)=L[n] and Pr(T)=1-L[n] then`):
print(`Alice and Bob have the same chances of winning. Try:`):
print(`LD({[1,1]},{[1,2]},100);`):

elif nops([args])=1 and op(1,[args])=Mihai then
print(`Mihai(m,A,B): inputs a positive integer m larger than 1, and two words in {1,..,m} of the same length A and B and outputs the number (according to Mihai Nica) k, such that`):
print(`asymptotically with n tosses of a fair coin, Pr(A>B)-Pr(A<B)  is asymptotic to k/sqrt(Pi)/sqrt(n). Try:`):
print(`Mihai(2,[1,1],[1,2]);`):

elif nops([args])=1 and op(1,[args])=MomsAB then
print(`MomsAB(m,A,B,n):  inputs m (the size of the alphabet, and sets A and B of words in {1,..,m} of the same lengths`):
print(`describing the Alice and Bob bets. Outpus the average occurrences`):
print(`The [[avA,VarA],[avB,VarB],co]: Try:`):
print(`MomsAB(2,{[1,1]},{[1,2]},n);`):


elif nops([args])=1 and op(1,[args])=NicaBook then
print(`NicaBook(m,n): inputs a positive integer m >=2 denoting that the alphabet is {1,..,m} and outputs the`):
print(`ranked lists of all counter-bets for all words in {1,...,m} with n letters, followed by the `):
print(`largest advantage. Try:`):
print(`NicaBook(2,3);`):

elif nops([args])=1 and op(1,[args])=NuSS then
print(`NuSS(w,Sv): Given a word w, and a set of words Sv, outputs the number of occurrence of members of Sv in w. Try:`):
print(`NuSS([1,1,1,2,2],{[1,1],[1,2]});`):


elif nops([args])=1 and op(1,[args])=NuSS1 then
print(`NuSS1(w,v): the number of consecutive subword (factor) v in the word w. For example, try:`):
print(`NuSS1([1,2,1,2,1,2,2],[1,2]);`):


elif nops([args])=1 and op(1,[args])=RECaz then
print(`RECaz(m,A,B,n,N): Given a positive integer`):
print(`m and two sets A,B, consisting  of words of the same length in the`):
print(`alphabet {1, ...,m}, `):
print(`uses the continuous Almkvist-Zeilberger algorithm,`):
print(`to find a linear recurrence operator with`):
print(`polynomial coefficients, ope(n,N) in terms of n and the`):
print(`shift operator in n, N, that annihilates`):
print(`the sequence`):
print(`"number of words in {1, ...m} with more occurrences of factors in A than factors in B"`):
print(`For example for the original Daniel Litt problem, type:`):
print(`RECaz(2,{[1,1]},{[1,2]},n,N);`):

elif nops([args])=1 and op(1,[args])=SeqFromRec then
print(`SeqFromRec(ope,Ini,n,N,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0, where ope(n,N) has order L,`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-2,[1],n,N,10);`):

elif nops([args])=1 and op(1,[args])=TestConjVar then
print(`TestConjVar(n,m,N): testing the conjecture that lower variance implies better chances,using n=N. Try:`):
print(`TestConjVar(2,3,50);`):


elif nops([args])=1 and op(1,[args])=TestMihai then
print(`TestMihai(m,n,K): does TestMihai1(m,A,B) for all pairs of words in {1,..,m} of length n that are not equivalent. Try:`):
print(`TestMihai(2,3,5000);`):

elif nops([args])=1 and op(1,[args])=TestMihai1 then
print(`TestMihai1(m,A,B,K): inputs a pos. integer m and two words A and B of the same length, compares the exact value of Mihai(m,A,B) with the estimate given by`):
print(`WhoWon(2,{A},{B},K); Try:`):
print(`TestMihai1(2,[1,1],[1,2],5000);`):



elif nops([args])=1 and op(1,[args])=TestNicaConj then
print(`TestNicaConj(m,n): tests Miahi Nica's conjecture about the limit of the third moment divide by the second moment for different words A and B of n letters in the alphabet {1,...,m} `):
print(`try:`):
print(`TestNicaConj(2,5):`):

elif nops([args])=1 and op(1,[args])=TestNicaConj1 then
print(`TestNicaConj1(m,A,B): tests Miahi Nica's conjecture about the limit of the third moment divide by the second moment for words A and B in the alphabet {1, ..., m} `):
print(`try:`):
print(`TestNicaConj1(2,[1,1],[1,2]):`):

elif nops([args])=1 and op(1,[args])=WhoWon then
print(`WhoWon(m,A,B,K): Given a positive integer m>=2, and two sets of words (all consisting of the same length) in the alphabet {1,...,m}, A,B, and a large positive integer K`):
print(`decides whether for all n, of you you toss a fair m-sided coin n times you are more likely to have more A's than B's or the other way round. It also outputs the`):
print(`estimates of the constants alpha, beta, such that the prob. that A wins is 1/2-alpha/sqrt(n) and the prob. that B wins is 1/2-beta/sqrt(n). Try:`):
print(`WhoWon(2,{[1,1]},{[1,2]},30000);`):


elif nops([args])=1 and op(1,[args])=WhoWonN then
print(`WhoWonN(m,A,B,K): Given a positive integer m>=2, and two sets of words (all consisting of the same length) in the alphabet {1,...,m}, A,B, and a positive integer K`):
print(`determines whether it is more likely to have more occurrences of the members of A than B, or vice versa, followed by the`):
print(`respective, numerical probabilities. For example for the original Daniel Litt problem with HT vs. HH and 100 coin tosses, type:`):
print(`WhoWonN(2,{[1,2]},{[1,1]},100);`):


elif nops([args])=1 and op(1,[args])=WhoWonV then
print(`WhoWonV(m,A,B,K,co): verbose form of WhoWon(m,A,B,K) (q.v.). It starts, "Question Number co". Try: `):
print(`WhoWonv(2,{[1,1]},{[1,2]},10000,1);`):


elif nops([args])=1 and op(1,[args])=WhoWonVar then
print(`WhoWonVar(m,A,B): Does A have a smaller variance than B? Try:`):
print(`WhoWonVar(2,{[1,1]},{[1,2]});`):

elif nops([args])=1 and op(1,[args])=Words then
print(`Words(n,m): the set of n-letter words in the alphabet {1,...,m}. Try:`):
print(`Words(5,2);`):

elif nops([args])=1 and op(1,[args])=WtE then
print(`WtE(n,m,A,B,a,b): The weight-enumerator of the set of n-letter words in {1, ...,m} according to the weight a^NuSS(w,A)*b^NuSS(w,B). Done by brute force. Only for checking.Try:`):
print(`WtE(7,2,{[1,1]},{[2,2]},a,b);`):

elif nops([args])=1 and op(1,[args])=WW then
print(`WW(m,A,B,K1,K2): Same as WhoWon(m,A,B,K) but much faster for moderate K2, and only returning true (if A won), false (if B won) and 0 (if they equally likely)`):
print(`by consistently looking at the scores from n=K1 and n=K2. If it is inconsistent, it returns FAIL. Try:`):
print(`WW(2,{[1,1]},{[1,2]},20,400);`):


else
print(`There is no ezra for`,args):
fi:
 
end:


#########################start from David_Ian
 

 
#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:

 
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:

#########################End from David_Ian


###Start from FindRec.txt
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


#if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
# RETURN(FAIL):
#fi:

#if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
# RETURN(FAIL):
#fi:

#if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
# RETURN(FAIL):
#fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:





 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:
###End from FindRec.txt

#Words(n,m): the set of n-letter words in the alphabet {1,...,m}. Try:
#Words(5,2);
Words:=proc(n,m) local gu,i,w:
option remember:

if n=0 then
 RETURN({[]}):
fi:

gu:=Words(n-1,m):

{seq(seq([op(w),i],i=1..m), w in gu)}:

end:

#NuSS1(w,v): the number of consecutive subword (factor) v in the word w. For example, try:
#NuSS1([1,2,1,2,1,2,2],[1,2]);

NuSS1:=proc(w,v) local i,c:
c:=0:

for i from 1 to nops(w)-nops(v)+1 do
 if [op(i..i+nops(v)-1,w)]=v then
   c:=c+1:
 fi:
od:
c:
end:

#NuSS(w,Sv): Given a word w, and a set of words Sv, outputs the number of occurrence of members of Sv in w. Try:
#NuSS([1,1,1,2,2],{[1,1],[1,2]});
NuSS:=proc(w,Sv) local v:
add(NuSS1(w,v), v in Sv):
end:


#AvsB(n,m,A,B): inputs pos. integers n and m, and set of words A and B in the alphabet {1,...,m} finds the number of n-letter words in {1,..,m} where
#there are more occurrences of members of A then B. Try:
#AvsB(5,2,{[1,1]},{[1,2]});
AvsB:=proc(n,m,A,B) local gu,w,c:
c:=0:
gu:=Words(n,m):

for w in gu do
 if NuSS(w,A)>NuSS(w,B) then
   c:=c+1:
 fi:
od:
c:
end:

#ABseqBF(n,m,A,B): the first n terms of the sequence (starting at n1=1) of the sequence AvsB(n1,m,A,B) by brute force. For checking only. Try:
#ABseqBF(10,2,{[1,1]},{[1,2]});
ABseqBF:=proc(n,m,A,B) local n1:
[seq(AvsB(n1,m,A,B) ,n1=1..n)]:
end:


#WtE(n,m,A,B,a,b): The weight-enumerator of the set of n-letter words in {1, ...,m} according to the weight a^NuSS(w,A)*b^NuSS(w,B). Try:
#WtE(7,2,{[1,1]},{[2,2]},a,b);
WtE:=proc(n,m,A,B,a,b) local gu,w:
gu:=Words(n,m):
add(a^NuSS(w,A)*b^NuSS(w,B), w in gu):
end:




#GFwtE(m,A,B,a,b,x): inputs a positive integer m, sets of words in {1,...,m}, A,B, variables, a and b, and another variable x, finds the rational function
#Sum(WtE(n,m,A,B,a,b)*x^n,n=0..infinity); Try:
#GFwtE(2,{[1,1]},{[1,2]},a,b,x);
GFwtE:=proc(m,A,B,a,b,x) local i,f,t,z:

f:=GJgfDetail({seq(i,i=1..m)},{op(A),op(B)},x,t):
f:=subs({seq(x[i]=x,i=1..m)},f):
f:=subs({seq(t[op(z)]=a,z in A)},f):
f:=subs({seq(t[op(z)]=b,z in B)},f):
normal(f):
end:


#GFwtEold(m,A,B,a,b,x): inputs a positive integer m, sets of words in {1,...,m}, A,B, variables, a and b, and another variable x, finds the rational function
#Sum(WtE(n,m,A,B,a,b)*x^n,n=0..infinity); Try:
#GFwtEold(2,{[1,1]},{[1,2]},a,b,x);
GFwtEold:=proc(m,A,B,a,b,x) local k,eq,var, A1,lu,eq1,i,ru,lu1,X,gu,n1:

option remember:

k:={seq(nops(A1), A1 in A), seq(nops(A1), A1 in B)}:
if nops(k)<>1 then
 print(`all words must be the same length`):
RETURN(FAIL):
fi:
k:=k[1]:

if k<2 then
 print(`The lengths of the words must be at least 2`):
 RETURN(FAIL):
fi:

lu:=Words(k-1,m):

var:={seq(X[lu1],lu1 in lu)}:

eq:={}:

for lu1 in lu do
eq1:=X[lu1]-a^NuSS(lu1,A)*b^NuSS(lu1,B)*x^nops(lu1):
for i from 1 to m do
  ru:=x:
 if member([op(lu1),i],A) then
  ru:=ru*a:
 fi:
 if member([op(lu1),i],B) then
  ru:=ru*b:
 fi:
eq1:=eq1 -ru*X[[op(2..nops(lu1),lu1),i]]:
od:
eq:=eq union {eq1}:
od:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

gu:=add(X[lu1],lu1 in lu):
gu:=subs(var,gu):


gu:=gu+ add(x^n1*WtE(n1,m,A,B,a,b),n1=0..k-2):
normal(gu):
end:


#CheckGFwtE(m,A,B,n): checks GFwtE(m,A,B,a,b,x) up to the n-th coefficient. Try:
#CheckGFwtE(2,{[1,1]},{[1,2]},10);
CheckGFwtE:=proc(m,A,B,n) local a,b,x,n1,f:
f:=GFwtE(m,A,B,a,b,x):

if f=FAIL then
 RETURN(FAIL):
fi:

f:=taylor(f,x=0,n+4):
evalb({seq(expand(coeff(f,x,n1)-WtE(n1,m,A,B,a,b)),n1=0..n)}={0}):
end:


#ABseq(n,m,A,B): the first n terms of the sequence (starting at n1=1) of the sequence AvsB(n1,m,A,B) done cleverly. Try:
#ABseq(10,2,{[1,1]},{[1,2]});
ABseq:=proc(n,m,A,B) local n1,x,a,b,f,t,L,f1,i1:
f:=GFwtE(m,A,B,a,b,x):

if f=FAIL then
 RETURN(FAIL):
fi:

f:=normal(subs({a=t,b=1/t},f)):

f:=taylor(f,x=0,n+10):

L:=[]:

for n1 from 1 to n do
f1:=coeff(f,x,n1):
L:=[op(L),add(coeff(f1,t,i1),i1=1..degree(f1,t))]:
od:

L:

end:


###Start from Almkvist-Zeilberger




####Start DiffToREc1


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:

#DiffToRec1(ope,t,Dt,n,N): Given a linear differential operator
#ope(t,Dt) and symbols n,N, outputs the linear recurrence
#operator with polynomial coefficients satisfied by
#the sequence of coeff. of any f.p.s. annihilating ope(t,Dt)
#where N is the shift in n. For example, try:
#DiffToRec1(Dt+t, t,Dt,n,N);
DiffToRec1:=proc(ope,t,Dt,n,N) local gu,i,j,ope1,d,j1:
ope1:=expand(ope):

gu:=0:

for i from ldegree(ope1,t) to degree(ope1,t) do
 for j from 0 to degree(ope1,Dt) do
   gu:=gu+coeff(coeff(ope1,t,i),Dt,j)*mul(n-i+j1,j1=1..j)*N^(j-i):
 od:
od:


d:=ldegree(gu,N):
gu:=N^(-d)*subs(n=n-d,gu):



Yafe(gu,N)[2]:

end:

 

####End DiffToRec1


#HSP(ope,n,N): the highest singularity of the recurrence
#operator ope(n,N). For example, try:
#HSP((n-3)*(n-4)*N+1,n,N);
HSP:=proc(ope,n,N) local ope1,pol,gu,i,gu1:
ope1:=numer(normal(ope)):
pol:=lcoeff(ope1,N):
gu:=[solve(pol,n)]:
gu1:={}:
for i from 1 to nops(gu) do
 if type(gu[i],integer) and gu[i]>=0 then
   gu1:=gu1 union {gu[i]}:
 fi:
od:

max(op(gu1)):
end:

#SeqFromRec(ope,Ini,n,N,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
#For example, try:
#SeqFromRec(N-2,[1],n,N,10);
SeqFromRec:=proc(ope,Ini,n,N,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):


if HSP(ope1,n,N)>=1 then
 RETURN(FAIL):
fi:

if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:
 
pashetc:=proc(p,D)
local gu,p1,i,mekh:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,D) do
gu:=gu+factor(coeff(p1,D,i))*D^i:
od:
RETURN(gu,mekh):
end:

gzor:=proc(f,x,n)
local i,gu:
gu:=f:
for i from 1 to n do
gu:=diff(gu,x):
od:
RETURN(gu):
end:


goremc:=proc(D,ORDER)
local i,gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*D^i:
od:
RETURN(gu):
end:

duisc:= proc(A,ORDER,y,x,D)
local KAMA,gorem,conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,i,shad:
 
KAMA:=50:
gorem:=goremc(D,ORDER):
 
conj:=gorem:
yakhas:=0:
 
for i from 0 to degree(conj,D) do
yakhas:=yakhas+normal(coeff(conj,D,i)*gzor(A,x,i)/A):
yakhas:=normal(yakhas):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve(eq,va1):
fu:=normal(subs(va1,fu)):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
shad:=pashetc(gorem,D):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
shad[1],S:
end:


Mygcd:=proc(L) local i,L1:
if nops(L)=1 then
 RETURN(1):
fi:
if nops(L)=2 then
 RETURN(gcd(L[1],L[2])):
fi:
L1:=[op(1..nops(L)-1,L)]:
gcd(Mygcd(L1),L[nops(L)]):
end:

AZc:=proc(A,y,x,D)
local ORDER,gu,KAMA,ope,C,kak,i1:
KAMA:=50:
 
for ORDER from 1 to KAMA do
gu:=duisc(A,ORDER,y,x,D):


if gu[1]<>0 then
ope:=gu[1]:
C:=gu[2]:
kak:=Mygcd([seq(coeff(ope,D,i1),i1=0..ORDER)]):
if kak<>1 then
ope:=add(normal(coeff(ope,D,i1)/kak)*D^i1,i1=0..ORDER):
C:=normal(C/kak):
fi:

 RETURN(ope,C):
fi:
od:
0:
end:

#RSct(m,w1,w2,z,s): The rational function in z and s
#such that the constant term of s is the generating
#function for words in {1,..,m}, according to length that
#have the same number of occurences of the factor w1
#and w2.
#For example, try:
#RSct(2,[1,1],[1,2],z,s);
RSct:=proc(m,w1,w2,z,s) local gu,t,i1:


if not type(m, integer) and m>=1 then
 print(`Bad input`):
 RETURN(FAIL):
fi:


if not (type(w1,list) and type(w2,list)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if nops(w1)<>nops(w2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if {op(w1)} minus {seq(i1,i1=1..m)}<>{} then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if {op(w2)} minus {seq(i1,i1=1..m)}<>{} then
 print(`Bad input`):
 RETURN(FAIL):
fi:


if not (type(w1,list) and type(w2,list) and 
{op(w1)} minus {seq(i1,i1=1..m)}={} and
{op(w2)} minus {seq(i1,i1=1..m)}={} and
w1<>w2) then
 print(`bad input`):
 RETURN(FAIL):
fi:

gu:=WtE(m,nops(w1),z,t,{w1,w2}):

normal(subs({t[op(w1)]=s,t[op(w2)]=1/s},gu)):

end:



#DOaz(m,A,B,z,Dz): Given a positive integer
#m and two sets of words A and B of words in {1,...,m} of the same length in the
#uses the continuous Almkvist-Zeilberger algorithm,
#to find a linear differential operator with
#polynomial coefficients, in terms of z and the
#differentiation w.r.t. z, Dz, that annihilates
#the (ordinary) generating function
#of the sequence
#"number of words in {1, ...m} of length n that have more occurrences of 
#words in A then words in B
#For example for the original Daniel Litt problem
#DOaz(2,{[1,1]},{[1,2]},t,Dt);
DOaz:=proc(m,A,B,z,Dz) local mu,gu,t,a,b,f:
option remember:

f:=GFwtE(m,A,B,a,b,z):
f:=normal(subs({a=1/t,b=t},f)/(1-t)):
if subs(t=0,f)=0 then
f:=t*f:
f:=normal(subs(t=1/t,f))/t
fi:

gu:=AZc(f,t,z,Dz):
if gu<>0 then
 RETURN(gu[1]):
else
 RETURN(FAIL):
fi:
end:

#RECaz(m,A,B,n,N): Given a positive integer
#m and two sets A,B, consisting  of words of the same length in the
#alphabet {1, ...,m}, 
#uses the continuous Almkvist-Zeilberger algorithm,
#to find a linear recurrence operator with
#polynomial coefficients, ope(n,N) in terms of n and the
#shift operator in n, N, that annihilates
#the sequence
#"number of words in {1, ...m} with more occurrences of factors in A than factors in B"
#For example for the original Daniel Litt problem, type:
#RECaz(2,{[1,1]},{[1,2]},n,N);
RECaz:=proc(m,A,B,n,N) local ope,z,Dz,ORDER,lu,i1:
option remember:
ope:=DOaz(m,A,B,z,Dz):
ope:=DiffToRec1(ope,z,Dz,n,N):
ORDER:=degree(ope,N):
ope:=add(factor(coeff(ope,N,i1))*N^i1,i1=0..ORDER):
lu:=ABseq(ORDER,m,A,B):

if ABseq(3*ORDER,m,A,B)<>SeqFromRec(ope,lu,n,N,3*ORDER) then
 print(`Something bad happened`):
  RETURN(FAIL):
fi:

[ope,lu]:
end:

###End from Almkvist-Zeilberger



#WhoWonN(m,A,B,K): Given a positive integer m>=2, and two sets of words (all consisting of the same length) in the alphabet {1,...,m}, A,B, and a positive integer K
#determines whether it is more likely to have more occurrences of the members of A than B, or vice versa, followed by the
#respective, numerical probabilities. For example for the original Daniel Litt problem with HT vs. HH and 100 coin tosses, type:
#WhoWonN(2,{[1,2]},{[1,1]},100);

WhoWonN:=proc(m,A,B,K) local a,b:
a:=evalf(ABseq(K,m,A,B)[K]/m^K):
b:=evalf(ABseq(K,m,B,A)[K]/m^K):
if a>b then
[true,a,b]:
elif a<b then
[false,a,b]:
else
[0,a,b]:
fi:
end:

#WhoWon(m,A,B,K): Given a positive integer m>=2, and two sets of words (all consisting of the same length) in the alphabet {1,...,m}, A,B, and a large positive integer K
#decides whether for all n, of you you toss a fair m-sided coin n times you are more likely to have more A's than B's or the other way round. It also outputs the
#estimates of the constants alpha, beta, such that the prob. that A wins is 1/2-alpha/sqrt(n) and the prob. that B wins is 1/2-beta/sqrt(n). Try:
#WhoWon(2,{[1,1]},{[1,2]},30000);
WhoWon:=proc(m,A,B,K) local n,N, ope1,ope2,L1,L2,alpha,beta,lu,d,i:
if K<1000 then
 print(`K has to be at least`, 1000):
 RETURN(FAIL):
fi:

if ABseq(100,m,A,B)=ABseq(100,m,B,A) then
  RETURN(0):
fi:


ope1:=RECaz(m,A,B,n,N):
ope2:=RECaz(m,B,A,n,N):
L1:=SeqFromRec(op(ope1),n,N,K):
L2:=SeqFromRec(op(ope2),n,N,K):

if L1=L2 then
 RETURN(0):
fi:

if min(op(10..K,L1-L2))>0 and max(op(10..K,L1-L2))>0 then
 lu:=true:
elif min(op(10..K,L2-L1))>0 and max(op(10..K,L2-L1))>0 then
 lu:=false:
else
 RETURN(FAIL):
fi:

L1:=evalf([seq((1/2-L1[i]/m^i)*sqrt(i),i=K-100..K)]):

d:=trunc(-log[10]( abs(L1[1]-L1[-1]))):
if d<3 then
 alpha:=FAIL:
else
alpha:=evalf(L1[-1],d-2):
fi:


L2:=evalf([seq((1/2-L2[i]/m^i)*sqrt(i),i=K-100..K)]):

d:=trunc(-log[10]( abs(L2[1]-L2[-1]))):
if d<3 then
 alpha:=FAIL:
else
beta:=evalf(L2[-1],d-2):
fi:

[lu,alpha,beta]:

end:


#WhoWonV(m,A,B,K,co): Verbose version of WhoWon(m,A,B,K). Try:
#WhoWonV(2,{[1,1]},{[1,2]},30000,1);
WhoWonV:=proc(m,A,B,K,co) local n,N, ope1,ope2,L1,L2,alpha,beta,lu,d,i,t0,ini1,ini2:
t0:=time():

if K<1000 then
 print(`K has to be at least`, 1000):
 RETURN(FAIL):
fi:
print(``):
print(`---------------------------------------------------------------------------------------------------`):
print(``):


print(`Question Mumber`, co, `: If Alice bets that there are more occurrences of consecutive substrings in`, A, `than in`, B, `and Bob, bets vice versa, who is more likely to win?`):

if ABseq(100,m,A,B)=ABseq(100,m,B,A) then
print(`Answer: They are equally likely to win`):
  RETURN(0):
fi:


ope1:=RECaz(m,A,B,n,N):
ope2:=RECaz(m,B,A,n,N):

ini1:=ope1[2]:
ini2:=ope2[2]:
ope1:=ope1[1]:
ope2:=ope2[1]:

 
if  WW(m,A,B,10,200) then
 print(`Answer: Alice is more likely to win.`):
else
 print(`Answer: Bob is more likely to win.`):
fi:


print(`Let's be more quantitative. We need the following two lemmas`):

print(` Lemma `,cat(co,`.`, 1), `: Let a(n) be the number of words of length n in the alphabet `, {seq(i,i=1..m)}, `with more occurences of`, A, `then `,B):
print(`a(n) satisfies the following linear recurrence equation with polynomial coefficients of order`, degree(ope1,N)):
print(``):
print(add(coeff(ope1,N,i)*a(n+i),i=0..degree(ope1,N))=0):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(a(i)=ini1[i],i=1..nops(ini1))):
print(``):

print(` Lemma `,  cat(co,`.`, 2), `: Let b(n) be the number of words of length n in the alphabet `, {seq(i,i=1..m)}, `with more occurences of`, B, `then `,A):
print(`b(n) satisfies the following linear recurrence equation with polynomial coefficients of order`, degree(ope2,N)):
print(``):
print(add(coeff(ope2,N,i)*b(n+i),i=0..degree(ope2,N))=0):
print(``):

print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(b(i)=ini2[i],i=1..nops(ini2))):
print(``):

print(`This enables us to compute the first`,K, `terms of these sequences`):

L1:=SeqFromRec(ope1,ini1,n,N,K):
L2:=SeqFromRec(ope2,ini2,n,N,K):


if min(op(10..K,L1-L2))>0 and max(op(10..K,L1-L2))>0 then
 lu:=true:
elif min(op(10..K,L2-L1))>0 and max(op(10..K,L2-L1))>0 then
 lu:=false:
else
 RETURN(FAIL):
fi:


print(``):
L1:=evalf([seq((1/2-L1[i]/m^i)*sqrt(i),i=K-100..K)]):

d:=trunc(-log[10]( abs(L1[1]-L1[-1]))):
if d<3 then
 alpha:=FAIL:
else
alpha:=evalf(L1[-1],d-2):
fi:


L2:=evalf([seq((1/2-L2[i]/m^i)*sqrt(i),i=K-100..K)]):

d:=trunc(-log[10]( abs(L2[1]-L2[-1]))):
if d<3 then
 alpha:=FAIL:
else
beta:=evalf(L2[-1],d-2):
fi:


print(`The probability of Alice winning the bet, as n grows larger is approximately`):
print(``):
print(1/2-alpha/sqrt(n)):
print(``):


print(`The probability of Bob winning the bet, as n grows larger is approximately`):
print(``):
print(1/2-beta/sqrt(n)):
print(``):

if lu=true then
print(`You can see that indeed Alice has a better chance, but not significantly so.`):
else
print(`You can see that indeed Bob has a better chance, but not significantly so.`):
fi:
print(``):
print(`-----------------------------------------`):
print(``):
print(`This took`, time()-t0, `seconds to generate`):
print(``):
[lu,alpha,beta]:
end:




#WW(m,A,B,K1,K2): Same as WhoWon(m,A,B,K) but much faster for moderate K, and only returning true (if A won), false (if B won) and 0 (if they equally likely)
#by consistently looking at the scores from n=K1 and n=K2. If it is inconsistent, it returns FAIL. Try:
#WW(2,{[1,1]},{[1,2]},20,400);
WW:=proc(m,A,B,K1,K2) local L1,L2:


L1:=ABseq(K2,m,A,B):

L2:=ABseq(K2,m,B,A):


if L1=L2 then
 RETURN(0):
fi:

if min(op(K1..K2,L1-L2))>0 and max(op(K1..K2,L1-L2))>0 then
 RETURN(true):
elif min(op(K1..K2,L2-L1))>0 and max(op(K1..K2,L2-L1))>0 then
 RETURN(false):
else
 RETURN(FAIL):
fi:

end:



#BETc(BET,m): All the equivalence classes of B of words of length nops(BET) for B under "having the same GFwtE(m,{BET},{B},a,b,z):
BETc:=proc(BET,m) local k,lu,gu,T,a,b,x,B,gu1:
k:=nops(BET):
lu:=Words(k,m) minus {BET}:

gu:={seq(GFwtE(m,{BET},{B},a,b,x), B in lu)}:


for gu1 in gu do
T[gu1]:={}:
od:

for B in lu do
T[GFwtE(m,{BET},{B},a,b,x)]:=T[GFwtE(m,{BET},{B},a,b,x)] union {B}:
od:

{seq(T[gu1],gu1 in gu)}:

end:

#Book1(m,BET,K): inputs a positive integer m (at least 2), and a positive integer k (at least 2) Bob's bet, BET, ( a word in the alphabet {1,...,m})
#and a large positive intetger K (at least 1000) outputs
#a book of the winning chances in the Litt game of all k-letter words in {1,...,m} versus BET. Try:
#Book1(2,[1,1,1,1],5000);
Book1:=proc(m,BET,K) local k,lu,i,t0,lu1,T,gu,ku,T1,ku1,i1,fu,n,a1,co:
k:=nops(BET):
t0:=time():
if not (type(m,integer) and type(k,integer) and type(K,integer) and m>=2 and k>=2 and K>=1000) then
print(`Bad input`):
RETURN(FAIL):
fi:

lu:=BETc(BET,m):
co:=0:
print(``):
print(`The Relative Probability of Winning a Daniel Litt style game when rolling a fair`, m, `sides die marked with`, {seq(i,i=1..m)}, `of number of occurrences of`, BET, `Versus the number of occurences`):
print(`of all`, k, `letter words in`,  {seq(i,i=1..m)}):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
ku:={}:

for lu1 in lu do
co:=co+1:
 gu:=WhoWonV(m,{lu1[1]},{BET},K,co):
if gu=0 then
 T[lu1]:=0:
ku:=ku union {0}:
else
 T[lu1]:=gu[3]-gu[2]:
ku:=ku union {T[lu1]}:
fi:
od:

for ku1 in ku do
T1[ku1]:={}:
od:
for lu1 in lu do
 T1[T[lu1]]:=T1[T[lu1]] union {lu1}:
od:


ku:=sort(convert(ku,list)):
fu:=[]:

for i1 from 1 to nops(ku) do
fu:=[op(fu),[T1[ku[i1]],ku[i1]]]:
od:

fu:=[seq(fu[nops(fu)+1-i1],i1=1..nops(fu))]:

if fu[1][2]<>0 then
print(`The best bet against`, BET, `is a member of`, {seq(op(a1),a1 in fu[1][1])}, `and then your edge, if you have n rolls, is approximately`, fu[1][2]/sqrt(n)):
else
print(`The best bet against`, BET, `is a member of`,  {seq(op(a1),a1 in fu[1][1])}, `and then your edge, if you have n rolls, is exactly 0`):
fi:
print(``):
for i from 2 to nops(fu) do
print(``):
if fu[i][2]<>0 then
 print(`The next best bet is a member of`,  {seq(op(a1),a1 in fu[i][1])}, `and then your edge is approximately`, fu[i][2]/sqrt(n)):
else
 print(`The next best bet is a member of`, {seq(op(a1),a1 in fu[i][1])}, `and then your edge is exactly 0`):
fi:
od:


print(`-----------------------`): 
print(``):
print(`This ends this chapter that took`, time()-t0, `seconds to generate. `):
print(``):

end:



#Temu(w,pi): the mapping of w under pi
Temu:=proc(w,pi) local i:
[seq(pi[w[i]],i=1..nops(w))]:
end:


#EWords(n,m): A set of representatives of  n-letter words in {1,...,m}
EWords:=proc(n,m) local mu,lu,gu,lu1:
mu:=Words(n,m):
lu:=permute(m):

gu:={}:

while mu<>{} do
 gu:=gu union {mu[1]}:
 mu:=mu minus {seq(Temu(mu[1],lu1),lu1 in lu)}:
od:
gu:
end:


#Book(m,k,K): A book about the relative probablity of winning Daniel Litt bets for all k-letter words in the alphabet {1,...,m} using K terms. Try:
#Book(2,3,5000);
Book:=proc(m,k,K) local gu,gu1,co,t0:
gu:=EWords(k,m):

t0:=0:
co:=0:
print(`How to bet against any given`, k, `-letter word in the Daniel Litt game with a fair die with`, m, ` faces `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`-----------------------------------------------------------------`):
print(``):

for gu1 in gu do
co:=co+1:
print(``):
print(`-----------------------------------------------------------------`):
print(``):
print(`Chapter Number`, co):
print(``):
Book1(m,gu1,K):
print(``):
print(``):
od:
print(`--------------------------------------------------------`):
print(``):
print(`This ends this book that took`, time()-t0, `seconds to generate. `):
print(``):
end:




#Info(m,A,B,r,K,n,N,C): Given a positive integer m>=2, and two sets of words (all consisting of the same length) in the alphabet {1,...,m}, A,B, and a large positive integer K and a positive integer
#r, outputs the pair of lists
#[[ope1,Ini1,Asy1],[ope2,Ini2,Asy2]]. 
#Where ope1 is the operator (of complexity at most C) annihilating the sequence 2^(n-1)-Number of n-letter words with more A's than B, Ini1 the initial condition and Asy1 the
#asympotics to order r, 
#[ope2,Ini2,Asy2] same but with 2^(n-1)-Number of n-letter words with more B's than A. Try:
#Info(2,{[1,1]},{[1,2]},5,5000,n,N,10);
Info:=proc(m,A,B,r,K,n,N,C) local K1, ope1,ope2,L1,L2,Ini1,Ini2,i,d1,o1,i1,ku1,ku2,c1,c2,as1,as2,c11,c21:

K1:=max(seq(seq((2+d1)*(1+o1)+10,o1=0..C-d1),d1=0..C))+40:

if ABseq(K1,m,A,B)=ABseq(K1,m,B,A) then
  RETURN(0):
fi:

L1:=ABseq(K1,m,A,B):
L2:=ABseq(K1,m,B,A):

L1:=[seq(L1[i]/m^i,i=1..nops(L1))]:
L2:=[seq(L2[i]/m^i,i=1..nops(L2))]:

L1:=[seq(1/2-L1[i],i=1..nops(L1))]:
L2:=[seq(1/2-L2[i],i=1..nops(L2))]:

ope1:=Findrec(L1,n,N,C):
if ope1=FAIL then
 RETURN(FAIL):
fi:

d1:=degree(numer(ope1),n):
o1:=degree(ope1,N):

for i1 from 1 to d1-1 while ope1<>FAIL do
ope1:=findrec(L1,d1-i1,o1+i1,n,N):
od:

ope2:=Findrec(L2,n,N,C):
if ope1=FAIL then
 RETURN(FAIL):
fi:

d1:=degree(numer(ope2),n):
o1:=degree(ope2,N):

for i1 from 1 to d1-1 while ope2<>FAIL do
ope2:=findrec(L2,d1-i1,o1+i1,n,N):
od:

Ini1:=[op(1..degree(ope1,N),L1)]:
Ini2:=[op(1..degree(ope2,N),L2)]:



read `AsyRecSilent.txt`:

ku1:=[AsyC(ope1,n,N,r,Ini1,K)]:

ku2:=[AsyC(ope2,n,N,r,Ini2,K)]:


c1:=ku1[1]:
c2:=ku2[1]:
as1:=expand(ku1[2]):
as2:=expand(ku2[2]):
c11:=simplify(op(1,as1)*sqrt(n),symbolic):

if not type(c11,numeric) then
 print(c11):
RETURN(FAIL):
fi:


c21:=simplify(op(1,as2)*sqrt(n),symbolic):


if not type(c21,numeric) then
 print(c21):
RETURN(FAIL):
fi:

c1:=c1*c11:
c2:=c2*c21:
as1:=expand(as1/c11):
as2:=expand(as2/c21):

if type(c1,numeric) then
 c1:=identify(evalf(c1*sqrt(Pi)))/sqrt(Pi):
fi:

if type(c2,numeric) then
 c2:=identify(evalf(c2*sqrt(Pi)))/sqrt(Pi):
fi:

RETURN([[ope1,Ini1,[c1,as1]], [ope2,Ini2,[c2,as2]]]):

end:


#MyFsolve(Eq,p)
MyFsolve:=proc(Eq,p) local gu,i:
gu:=[fsolve(Eq,p)]:

if gu=[] then
 RETURN(FAIL):
fi:

for i from 1 to nops(gu) do
if gu[i]>0 and gu[i]<1 then
  RETURN(gu[i]):
fi:
od:
FAIL:
end:



#LD(A,B,N): inputs two sets of words A and B consisting of words all of the same
#length, and a positive integer N outputs the list of 
#length N such that if Pr(H)=L[n] and Pr(T)=1-L[n] then
#Alice and Bob have the same chances of winning. Try:
#LD({[1,1]},{[1,2]},100);
LD:=proc(A,B,N) local p,f,x,t,a,b,f1A,f1B,L,n,f1,i,a1,ka:

#if {op(A),op(B)}<>{1,2} then
# print(A, B, `must be subsets of`, {1,2}):
# RETURN(FAIL):
#fi:

if nops({seq(nops(a1), a1 in A),seq(nops(a1),a1 in B)})<>1 then
print(`all words must be the same length`):
RETURN(FAIL):
fi:

f:=GJgfDetail({1,2},{op(A),op(B)},x,t);


f:=subs({seq(t[op(a1)]=a,a1 in A),seq(t[op(a1)]=b,a1 in B)},f):


f:=subs({x[1]=p*x,x[2]=(1-p)*x},f);
f:=subs({a=t,b=1/t},f):
f:=taylor(f,x=0,N+3):
L:=[0$10]:
for n from 11 to N do
 f1:=coeff(f,x,n):
 f1A:=expand(add(coeff(f1,t,i),i=1..degree(f1,t))):
 f1B:=expand(add(coeff(f1,t,i),i=ldegree(f1,t)..-1)):
 ka:=MyFsolve(f1A=f1B,p):
 L:=[op(L),ka]:
od:
 

L:
end:



#CoeffXn(f,x,n): given a rational function whose denominator is a power of (1-x) finds the coeff. of x^n in its Maclaurin expansion
CoeffXn:=proc(f,x,n) local f1,r,t,i:
r:=rem(numer(f),denom(f),x):
f1:=r/denom(f):
f1:=normal(subs(x=1-1/t,f1)):
if degree(denom(f1),t)<>0 then
  RETURN(FAIL):
fi:
expand(add(coeff(f1,t,i)*binomial(n+i-1,i-1),i=1..degree(f1,t))):
end:


#MomsAB(m,A,B):  inputs m (the size of the alphabet, and sets A and B of words in {1,..,m} of the same lengths
#describing the Alice and Bob bets. Outpus the average occurrences
#The [[avA,VarA],[avB,VarB],co]: Try:
#MomsAB(2,{[1,1]},{[1,2]},x);
MomsAB:=proc(m,A,B,n) local a,b,f,av1,av2,m1,m2,co,x:
f:=GFwtE(m,A,B,a,b,x):
f:=subs(x=x/2,f):
av1:=CoeffXn(normal(subs({a=1,b=1},diff(f,a))),x,n):
av2:=CoeffXn(normal(subs({a=1,b=1},diff(f,b))),x,n):
m1:=expand(CoeffXn(normal(subs({a=1,b=1},a*diff(a*diff(f,a),a))),x,n)-av1^2):
m2:=expand(CoeffXn(normal(subs({a=1,b=1},b*diff(b*diff(f,b),b))),x,n)-av2^2):
co:=expand(CoeffXn(normal(subs({a=1,b=1},b*diff(a*diff(f,a),b))),x,n)-av1*av2):
[[av1,m1],[av2,m2],co]:
end:



#WhoWonVar(m,A,B): Does A have a smaller variance than B? Try:
#WhoWonVar(2,{[1,1]},{[1,2]});
WhoWonVar:=proc(m,A,B) local gu,n,m1,m2:
gu:=MomsAB(m,A,B,n):
m1:=gu[1][2]:
m2:=gu[2][2]:

evalb(lcoeff(m1,n)<lcoeff(m2,n)):

end:



#TestConjVar(n,m,N): testing the conjecture that lower variance implies better chances,using n=N. Try:
#TestConjVar(2,3,50);
TestConjVar:=proc(n,m,N) local gu,i,j:
gu:=convert(Words(n,m),list):

for i from 1 to nops(gu) do
 for j from i+1 to nops(gu) do
  if WhoWonVar(m,{gu[1]},{gu[2]})<>WhoWonN(m,{gu[1]},{gu[2]},N)[1] then
   print([gu[i],gu[j]]):
   RETURN(false):
  fi:
od:
od:
true:
end:



#GP1(L,n,d): Inputs a list of numbers L and a variable n and a non-neg. d outputs the polynomial of degree d
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP1:=proc(L,n,d) local a,i,P,eq,var:
if nops(L)<d+2 then
 ERROR(`List too short`):
 fi:
#P=a[0]+a[1]*n+...+a[d]*n^d
P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(L[i]=subs(n=i,P),i=1..d+2)}:
var:=solve(eq,var):
 if var=NULL then
   RETURN(FAIL):
 fi:
P:=subs(var,P):
if {seq(L[i]-subs(n=i,P),i=d+3..nops(L))}<>{0} then
  RETURN(FAIL):
 fi:
P:
end:



#GP(L,n): Inputs a list of numbers L and a variable n and outputs the polynomial of degree <=nops(L)-2
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP:=proc(L,n) local d,P:
for d from 0 to nops(L)-2 do
 P:=GP1(L,n,d):
 if P<>FAIL then
    RETURN(P):
 fi:
od:
FAIL:
end:

#Mihai(m,A,B): inputs two words in {0,1} of the same length A and B and outputs the number (according to Mihai Nica) k, such that
#asymptotically with n tosses of a fair coin, Pr(A>B)-Pr(A<B) =k/sqrt(Pi)/sqrt(n). Try:
#Mihai(2,[1,1],[1,2]);
Mihai:=proc(m,A,B) local n,f,x,a,b,t,m2,f1,f2,i,f3,m3:

if not (nops(A)=nops(B) and {op(A),op(B)}={seq(i,i=1..m)}) then
 RETURN(FAIL):
fi:

f:=GFwtE(m,{A},{B},a,b,x):
f:=normal(subs({a=t,b=1/t,x=x/m},f)):
f1:=t*diff(f,t):
if normal(subs(t=1,f1))<>0 then
print(normal(subs(t=1,f1))):
 print(`Something bad happened`):
 RETURN(FAIL):
fi:

f2:=normal(t*diff(f1,t)):

m2:=normal(subs(t=1,f2)):

m2:=GP([seq(coeff(taylor(m2,x=0,44),x,i),i=20..30)],n):

if m2=FAIL then
RETURN(FAIL):
fi:

m2:=normal(subs(n=n-19,m2)):

f3:=normal(t*diff(f2,t)):
m3:=normal(subs(t=1,f3)):

m3:=GP([seq(coeff(taylor(m3,x=0,44),x,i),i=20..30)],n):

if m3=FAIL then
 RETURN(FAIL):
fi:

m3:=normal(subs(n=n-19,m3)):

#-1/3*sqrt(limit(m3^2*n/m2^3,n=infinity))/sqrt(2):

-1/3*limit(m3*sqrt(n)/m2^(3/2),n=infinity)/sqrt(2):
end:




#TestMihai1(m,A,B,K): inputs a pos. integer m and two words A and B of the same length, compares the exact value of Mihai(m,A,B) with the estimate given by
#WhoWon(2,{A},{B},K); Try:
#TestMihai1(2,[1,1],[1,2],5000);
TestMihai1:=proc(m,A,B,K) local lu1,lu2:
 lu1:=evalf(Mihai(m,A,B)/sqrt(Pi)):
 lu2:=WhoWon(m,{A},{B},K):
[lu1,evalf(lu2[3]-lu2[2])]:
end:



#TestMihai(m,n,K): does TestMihai1(m,A,B) for all pairs of words in {1,..,m} of length n that are not equivalent. Try:
#TestMihai(2,3,5000);
TestMihai:=proc(m,n,K) local gu,i,j:
gu:=convert(Words(n,m),list):

for i from 1 to nops(gu) do
 for j from i+1 to nops(gu) do
 if WhoWonN(m,{gu[i]},{gu[j]},50)[1]<>0 and Mihai(m,gu[i],gu[j])<>FAIL then
  print(`doing `, gu[i], gu[j] ):
  lprint(TestMihai1(m,gu[i],gu[j],K)):
 fi:
 od:
od:
end:

#CounterBet(m,A): inputs a word A in {1,..,m} and outputs the randked list  of the counterbets, from best to worst. Try:
#CounterBets(2,[1,1,1]);
CounterBets:=proc(m,A) local gu,n,B,R,ka1,T,i:
n:=nops(A):
gu:=Words(n,m) minus {A}:

R:=[]:
for i from 1 to nops(gu) do
 B:=gu[i]:
 if WhoWonN(m,{A},{B},50)[1]<>0 and Mihai(m,A,B)<>FAIL then
   ka1:= evalf(Mihai(m,A,B)):
  if not member(ka1,{op(R)}) then
  R:=[op(R),ka1]:
    T[ka1]:={B}:

  else
  T[ka1]:=T[ka1] union {B}:
 fi:
fi:
od:
R:=sort(R,`>`):
[seq([T[R[i]],R[i]],i=1..nops(R))]:
end:


 
  
#NicaBook(m,n): inputs a positive integer m >=2 denoting that the alphabet is {1,..,m} and outputs the
#ranked lists of all counter-bets for all words in {1,...,m} with n letters, followed by the 
#largest advantage. Try:
#NicaBook(2,3);
NicaBook:=proc(m,n) local t0,gu,aluf, si,i,lu,hal,N:

t0:=time():

print(`Using Mihai Nica's Asymptotic formula to find  Bob's worst counterbets (hence best for Alice) for any Alice word in`, {seq(i,i=1..m)}, ` with `, n , ` letters  `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

gu:=convert(Words(n,m),list):

si:=0:

for  i from 1 to nops(gu) do
print(``):
print(`For the word`, gu[i]):
print(``):
print(`The ranked list of counterbets, and their values are`):
print(``):
lu:=CounterBets(m,gu[i]):
lprint(lu):
hal:=lu[1][2]:
if hal>si then
 si:=hal:
aluf:=[gu[i],lu[1][1]]:
fi:
od:

print(``):

print(`The most unfair bet is if Alice bets`, aluf[1], ` and Bob bets any word in the set`, aluf[2]):
print(``):
lprint(`The advantage then is`, evalf(si), `that means that Pr(A>B)-Pr(B>A) with N tosees is asymptotic to`, evalf(si/sqrt(Pi))/sqrt(N) ):
print(``):
print(``):
print(`----------------------`):
print(``):
print(`This took`, time()-t0, `seconds. `):
end:



#HO(m,w): the Hanging Overlap of the word w in an alpahbet of size m, according to Mihai Nica. Try:
#HO(2,[1,2,1,2,1]);
HO:=proc(m,w) local i,c,n:
n:=nops(w):
c:=0:
for i from 2 to n do
 if [op(i..n,w)]=[op(1..n-i+1,w)] then
  c:=c+1/m^(i-1):
 fi:
od:
c:
end:

#TestNicaConj1(m,A,B): tests Miahi Nica's conjecture about the limit of the third moment divide by the second moment for words A and B
#try:
#TestNicaConj1(2,[1,1],[1,2]):
TestNicaConj1:=proc(m,A,B) local n,f,x,a,b,t,m2,f1,f2,i,f3,m3:

if not (nops(A)=nops(B) and {op(A),op(B)}={seq(i,i=1..m)}) then
 RETURN(FAIL):
fi:

f:=GFwtE(m,{A},{B},a,b,x):
f:=normal(subs({a=t,b=1/t,x=x/m},f)):
f1:=t*diff(f,t):
if normal(subs(t=1,f1))<>0 then
print(normal(subs(t=1,f1))):
 print(`Something bad happened`):
 RETURN(FAIL):
fi:

f2:=normal(t*diff(f1,t)):

m2:=normal(subs(t=1,f2)):

m2:=GP([seq(coeff(taylor(m2,x=0,44),x,i),i=20..30)],n):

if m2=FAIL then
RETURN(FAIL):
fi:

m2:=normal(subs(n=n-19,m2)):

f3:=normal(t*diff(f2,t)):
m3:=normal(subs(t=1,f3)):

m3:=GP([seq(coeff(taylor(m3,x=0,44),x,i),i=20..30)],n):

if m3=FAIL then
 RETURN(FAIL):
fi:

m3:=normal(subs(n=n-19,m3)):

evalb(limit(m3/m2,n=infinity)=3*(HO(m,A)-HO(m,B))):


end:

#TestNicaConj(m,n): does TestNicaConj1(m,A,B) for all pairs of words in {1,..,m} of length n that are now equivalent. Try:
#TestNicaConj(2,5);
TestNicaConj:=proc(m,n) local gu,i,j:
gu:=convert(Words(n,m),list):

for i from 1 to nops(gu) do
 for j from i+1 to nops(gu) do
 if not TestNicaConj1(m,gu[i],gu[j]) then
  print([gu[i],gu[j]], `is a counterexample`):
  RETURN(false):
 fi:
 od:
od:
true:

end: