######################################################################
##AMP.txt: Save this file as  AMP.txt                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read AMP.txt<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################

 
#Created: Jan. 2020
 
print(`Created:  Jan. 2020`):
print(` This is AMP.txt `):
print(`It is one of the packages that accompany the article `):
print(`Exploring the Absent-Minded Passengers Problem`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

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

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



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

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(`---------------------------------------`):
 print(`For a list of the Findrec procedures type ezraFindrec();, for help with`):
 print(`a specific procedure, type ezraFindrec(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):


ezraS:=proc()

if args=NULL then
 print(` The STORY procedures are:  Paper1, Paper2, Paper3 `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Anw1, Anw1Fa, Childn, ChildnWW, Cyc, FC1, Fstate, FstateG, GuessPol, GuessPol1G, Hnm, IniState, Pn, Wrn, Wt1, Wtpi `):
 print(``):

else
ezra(args):
fi:

end:

ezra1G:=proc()

if args=NULL then
 print(` The supporting generalized procedures are:  `):
 print(`  Anw1FoG,  Anw2FoG, Anw1G, Anw1FaG, ChildnWWg ,CoeffS, DerG, erk, PeelG, PeelG1, WtpiG,  Wt1G, WrnG,`):
 print(``):

else
ezra(args):
fi:

end:

ezraG:=proc()

if args=NULL then
 print(` The general weight procedures are: `):
 print(` AnwkG, AnwkGE, AnwkGsim  `):
 print(``):

else
ezra(args):
fi:

end:




ezra:=proc()

if args=NULL then
 print(`The main procedures are: AllStates,  AlphaPask1Asy, Anw1Fo ,   Anwk, AnwkE, AnwkSeq, AnwkSeqF, AnwkSim,  CoeffS, EVLL, MomPas, MomPask1, `):
  print(`MomPask1Asy, Opersk, Simk, SimkV, PrS, Vark  `):
 print(` `):



elif nops([args])=1 and op(1,[args])=AllStates then
print(`AllStates(n,k): The set of all states coming from n passengers {1, ....,n} the first k of them being absent-minded. Try:`):
print(`AllStates(5,1);`):

elif nops([args])=1 and op(1,[args])=AlphaPask1Asy then
print(`AlphaPask1Asy(n,r,K): The asymptotic expression for  the scaled r-th moment about the  mean of W_{n,1} up to K terms`):
print(`if r is even, and its square it r is odd.`):
print(`AlphaPask1Asy(n,4,6);`):

elif nops([args])=1 and op(1,[args])=Anw1 then
print(`Anw1(n,w): The full probability generating function for n passengers, the first being absent-minded,from the definition, By brute force`):
print(`Warning: very slow, for checking purposes only`):
print(`summing the weights of all possible permutations. Try:`):
print(`[seq(Anw1(n1,w),n1=1..10)]; `):



elif nops([args])=1 and op(1,[args])=Anw1Fa then
print(`Anw1Fa(n,w): The full probability generating function for n passengers, the first being absent-minded, using dynamical programming.`):
print(`much faster than Anw1(n,w). Try:`):
print(`[seq(Anw1Fa(n1,w),n1=1..10)]; `):

elif nops([args])=1 and op(1,[args])=Anw1Fo then
print(`Anw1Fo(n,w): The full probability generating function for n passengers, the first being absent-minded, using the explicit formula.`):
print(` Even faster  than Anw1Fa(n,w). Try:`):
print(`[seq(Anw1Fo(n1,w),n1=1..10)]; `):


elif nops([args])=1 and op(1,[args])=Anw1G then
print(`Anw1G(n,w): The full probability generating function for n passengers, the first being absent-minded,from the definition, By brute force`):
print(`Warning: very slow, for checking purposes only`):
print(`summing the GENERALIZED weights of all possible permutations.  The weight of a permutation is the product of w[i] over all i`):
print(`that are not in the i-th location, times the probability that it happens. `):
print(`Try:`):
print(`[seq(Anw1G(n1,w),n1=1..10)]; `):

elif nops([args])=1 and op(1,[args])=Anw1FaG then
print(`Anw1FaG(n,w): The full probability generating function for n passengers, the first being absent-minded,from the definition, By Dynamical programming.`):
print(`much faster than Anw1G(n,w)`):
print(`summing the GENERALIZED weights of all possible permutations.  The weight of a permutation is the product of w[i] over all i`):
print(`that are not in the i-th location, times the probability that it happens. `):
print(`Try:`):
print(`[seq(Anw1FaG(n1,w),n1=1..10)]; `):

elif nops([args])=1 and op(1,[args])=Anw1FoG then
print(`Anw1FoG(n,w): The full probability generating function for n passengers, the first being absent-minded, By the`):
print(`explicit formula.`):
print(`even faster than Anw1FaG(n,w)`):
print(`summing the GENERALIZED weights of all possible permutations.  The weight of a permutation is the product of w[i] over all i`):
print(`that are not in the i-th location, times the probability that it happens. `):
print(`Try:`):
print(`[seq(Anw1FoG(n1,w),n1=1..10)]; `):


elif nops([args])=1 and op(1,[args])=Anw2FoG then
print(`Anw2FoG(n,w): The full probability generating function for n passengers, the FIRST and SECOND being absent-minded,from the definition,`):
print(` by the explicit formula.`):
print(`Try:`):
print(`[seq(Anw2FoG(n1,w),n1=1..10)]; `):


elif nops([args])=1 and op(1,[args])=Anwk then
print(`Anwk(n,w,k): the generating function in w whose coefficient of w^i is the probability that exactly i passengers will seat in `):
print(`a different seat that has been assigned to them of there are n passengers {1, ..., n} the first k of whom are absent-minded.`):
print(`USING DYNAMILCAL PROGRAMMING`):
print(`Try: Anwk(5,w,1);`):


elif nops([args])=1 and op(1,[args])=AnwkE then
print(`AnwkE(n,w,k): the generating function in w whose coefficient of w^i is the probability that exactly i passengers will seat in `):
print(`a different seat that has been assigned to them of there are n passengers {1, ..., n} the first k of whom are absent-minded.`):
print(`USING THE EXPLICIT FORMULA`):
print(`Try: AnwkE(5,w,1);`):




elif nops([args])=1 and op(1,[args])=AnwkG then
print(`AnwkG(n,w,k): the generating function in w[1], ..., w[n] whose coefficient of w[i1]..w[ir] is the probability that  `):
print(`passengers i1, ..., ir are seating in a different seat that has been assigned to them`):
print(`if there are n passengers {1, ..., n} the first k of whom, namely {1, ...k} are absent-minded. Try: `):
print(`AnwkG(5,w,1);`):


elif nops([args])=1 and op(1,[args])=AnwkGE then
print(`AnwkGE(n,w,k): Same as AnwkG(n,w,k) using the EXPLICIT formula. Try`):
print(` AnwkGE(5,w,1);`):

elif nops([args])=1 and op(1,[args])=AnwkGsim then
print(`AnwkGsim(n,w,k,N): simulates the absent-minded passenger with k absent minded passengers and a total of n passengers `):
print(`detailed weight`):
print(`where the first k ones are absent-minded`):
print(`N times and outputs`):
print(`the fraction of those where pi[n]<>n and the generating function according to the location of the non-fixed points`):

elif nops([args])=1 and op(1,[args])=AnwkSeq then
print(`AnwkSeq(N,w,k): The first N terms of the sequence Anwk(n,w,k). Try:`):
print(`AnwkSeq(40,w,2);`):


elif nops([args])=1 and op(1,[args])=AnwkSeqF then
print(`AnwkSeqF(N,w,k): The first N terms of the sequence Anwk(n,w,k). FAST VERSION for 1<=k<=4. Try:`):
print(`AnwkSeqF(40,w,2);`):

elif nops([args])=1 and op(1,[args])=AnwkSim then
print(`AnwkSim(n,w,k,N): simulates the absent-minded passenger with k absent minded passengers and a total of n passengers `):
print(`simple weight`):
print(`where the first k ones are absent-minded`):
print(`N times and outputs`):
print(`the fraction of those where pi[n]<>n and the generating function according to the location of the bnon-fixed points`):


elif nops([args])=1 and op(1,[args])=Childn then
print(`Childn(St): Given a state St=[Pas,Seats,AM] where Pa is the current list of passengers (in order), Seats is the current list of seats (in order)`):
print(`and AM is the subset of  Pas that is abent-minded, outputs the set of children, i.e. the set of possible states`):
print(`that results from seating the first passender in Pa. Try:`):
print(`Childn([{1,2,3,4},{1,2,3,4},{1}]);`):

elif nops([args])=1 and op(1,[args])=ChildnWW then
print(`ChildnWW(St,w): Given a state St=[Pas,Seats,AM] where Pa is the current list of passengers (in order), Seats is the current list of seats (in order)`):
print(`and AM is the subset of  Pas that is abent-minded, and a variable name,`):
print(`outputs the set of pairs [child, weight], where weight is w to the power [number of passengers seating not in their seats], for the set of possible states`):
print(`that results from seating the first passender in Pa. Try:`):
print(`ChildnWW([{1,2,3,4},{1,2,3,4},{1}],w);`):

elif nops([args])=1 and op(1,[args])=ChildnWWg then
print(`ChildnWWg(St,w): Given a state St=[Pas,Seats,AM] where Pa is the current list of passengers (in order), Seats is the current list of seats (in order)`):
print(`and AM is the subset of  Pas that is abent-minded, and a variable name,`):
print(`outputs the set of pairs [child, weight], where weight is the generalized weight with w[1], ... w[n], for the set of possible states`):
print(`that results from seating the first passender in Pa. Try:`):
print(`ChildnWWg([{1,2,3,4},{1,2,3,4},{1}],w);`):


elif nops([args])=1 and op(1,[args])=CoeffS then
print(`CoeffS(n,k,S): Inputs positive integers n and k, and a subset S={i1, ..., ik}, of {1,...n}, outputs the coefficient of`):
print(`w[i1]...w[ik] in AnwkG(n,w,k); Try:`):
print(`CoeffS(10,2,{3,10});`):

elif nops([args])=1 and op(1,[args])=Cyc then
print(`Cyc(pi): the list of cycles of the permutation pi, except for fixed points, try:`):
print(`Cyc([3,1,2,4]);`):

elif nops([args])=1 and op(1,[args])=DerG then
print(`DerG(F,w,S): inputs a polynomial F in w[1], ..., w[n] (with n>=k, does not matter, n it is not part of the input) and`):
print(`and a subset S of {1,...k}, outputs subs({w[i]=1, in in S},F) and`):
print(`F-mul(w[i], in in S)* subs({w[i]=1, i in S)},F). Try:`):
print(`DerG( AnwkG(5,w,1),w,{1,2});`):

elif nops([args])=1 and op(1,[args])=erk then
print(`erk(r,k,w): e_r(w_1,...,w_k) the weight enumerator of subsets of {1,...k} with r elements where the weight is`):
print(`mul(1-w[i],i in S)*mul(w[i],i in not S). Try:`):
print(`erk(2,4,w);`):

elif nops([args])=1 and op(1,[args])=EVLL then
print(`EVLL(n,k): the expectation and variance according to Norbert Henze and Guenter Last. Try:`):
print(`EVLL(10,2); `):

elif nops([args])=1 and op(1,[args])=FC1 then
print(`FC1(pi,i): the cycle belonging to i in the permutation pi. Try: `):
print(`FC1([2,3,1],1);`):


elif nops([args])=1 and op(1,[args])=Fstate then
print(`Fstate(St,w): the weight of  a state St. Try:`):
print(`Fstate([{1,2,3,4},{1,2,3,4},{1}],w);`):

elif nops([args])=1 and op(1,[args])=FstateG then
print(`FstateG(St,w): the generalized weight of  a state St. Try:`):
print(`FstateG([{1,2,3,4},{1,2,3,4},{1}],w);`):

elif nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(LL,L1,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and`):
print(`a positive integer r, outputs a polynomial P such that `):
print(`L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Try:`):
print(` GuessPol([[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]], ApplyPolSeq(x[1]+x[2],x,[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]),x);`):


elif nops([args])=1 and op(1,[args])=GuessPol1G then
print(`GuessPol1G(LL,L1,A,r,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and`):
print(`a positive integers A and r, outputs a polynomial P of degree<=A in x[1] and total degree <=r in {x[2], ..., x[r]}, if it exists, such that `):
print(`L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Otherwise it returns FAIL. Try:`):
print(` GuessPol1G([[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]], ApplyPolSeq(x[1]^3+x[2],x,[[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]]),3,1,x);`):

elif nops([args])=1 and op(1,[args])=Hnm then
print(`Hnm(n,m):=add(1/i^m,i=1..n-1). Try:`):
print(`Hnm(10,2);`):

elif nops([args])=1 and op(1,[args])=IniState then
print(`IniState(n,k): the initial state with n passengers and k absent-minded ones. Try:`):
print(`IniState(6,2);`):

elif nops([args])=1 and op(1,[args])=MomPas then
print(`MomPas(k,n,Hn,r): Finds an explicit expression for the expected number (if r=1)`):
print(`and for the r-th moment about the mean for r>=2 of displaced passengers with k absent-minded passengers. Try:,`):
print(`in terms of n and the Harmonic numbers Hn. Try:`):
print(`MomPas(1,n,Hn,1);`):

elif nops([args])=1 and op(1,[args])=MomPask1 then
print(`MomPask1(n,Hn,r): Guesses an expression for the r-th moment about the  mean of W_{n,1}. Try:`):
print(`MomPask1(n,Hn,4);`):

elif nops([args])=1 and op(1,[args])=MomPask1Asy then
print(`MomPask1Asy(n,r,K): The asymptotic expression for  the r-th moment about the  mean of W_{n,1} up to the K-th order. Try:`):
print(`MomPask1Asy(n,4,6);`):

elif nops([args])=1 and op(1,[args])=Opers then
print(`Opers(k,w,n,N): The linear recurrence equation annihilating the sequence Anwk(n,w,k) for k=1,2,3, starting at n=k. Try:`):
print(`Opers(3,w,n,N)`):


elif nops([args])=1 and op(1,[args])=Paper1 then
print(`Paper1(w,n): inputs symbols w and n and outputs an article with recurrences for`):
print(` the probability generating function for F_k(w,n+k-1): the probability generating function`):
print(`according to the number of passengers sitting in a wrong seat if there are n+k-1 passengers, k of whom are absent-minded, for`):
print(`k from 1 to 4. Try:`):
print(`Paper1(w,n):`):

elif nops([args])=1 and op(1,[args])=Paper2 then
print(`Paper2(n,Hn,M1,M2): inputs symbols n, H, and positive integers M1, and M2, outputs explicit expressions for`):
print(`the r-th moment about the mean, for r between 1 and M1, of the random variable "Number of Passengers seating in the wrong seat" where`):
print(`the FIRST passenger is absent-minded, as well as asymptotic expressions up to the M2-th order. Here`):
print(`Hn[r] is short for sum(1/i^r,i=1..n-1). Try:`):
print(`Paper2(n,Hn,4,6);`):

elif nops([args])=1 and op(1,[args])=Paper3 then
print(`Paper3(k,n,N,R): compares the statistical data up to the M1-th moment the result of simulating the Absent-Minded Passenger`):
print(`problem with the first k passengers being absent-minded with the actual thing using the exact probablity generating function`):
print(`Anwk(N1,w,k). Try:`):
print(`Paper3(4,30,1000,4);`):

elif nops([args])=1 and op(1,[args])=PeelG then
print(`PeelG(w,n,k): Finds the coefficients of mul(w[i],i in S)*mul(1-w[i],i not in S) for all subsets S of {1, ..., k}  of`):
print(`AnwkG(w,n,k).  Try:`):
print(`PeelG(w,7,2);`):

elif nops([args])=1 and op(1,[args])=PeelG1 then
print(`PeelG1(w,n,k,SeqS): Inputs a symbols w, and positive integers n and a sequence of subsets of {1,...k}, outputs what happens to`):
print(`AnwkG(n,w,k) when you keep "peeling it", by applying DerG(F,w,SeqS[i]) to the list SeqS. Try: `):
print(` PeelG1(w,7,2,[{1,2},{2},{1}]); `):


elif nops([args])=1 and op(1,[args])=Pn then
print(`Pn(n): The set of lists [1,i1, .., ir] such that 1<=i1<i2<...<ir<=n. Corresponding to the`):
print(`scenarios in the case where only the first passenger  is absent-minded. Try: `):
print(`Pn(5); `):
print(``):


elif nops([args])=1 and op(1,[args])=PrS then
print(`PrS(n,k,S): Inputs positive integers n and k, and a subset S of {1,...n}, outputs the probability`):
print(`that with n passengers, the first k of whom are absent-minded, the members of S will all sit in the wrong seat. Try:`):
print(`PrS(10,1,{10});`):

elif nops([args])=1 and op(1,[args])=Simk then
print(`Simk(k,n): Outputs a random permutation obtained from the absent-minded passenger problem with n passengers, the first k of whom`):
print(`are absent-minded. Try:`):
print(`Simk(3,20);`):

elif nops([args])=1 and op(1,[args])=SimkV then
print(`Simk(k,n): Outputs a random permutation obtained from the absent-minded passenger problem with n passengers, the first k of whom`):
print(`are absent-minded. Verbose version, also prints out intermediate steps. Try:`):
print(`SimkV(3,20);`):


elif nops([args])=1 and op(1,[args])=Vark then
print(`Vark(k,n,Hn): The explicit expression for the variance according to Hence and last for W_{n,k}. Try: `):
print(`Vark(2,n,Hn)`):

elif nops([args])=1 and op(1,[args])=Wrn then
print(`Wrn(r,n,w): the weight enumerator of all sequences L that end in r. Using dynamical programing. Try:`):
print(`Wrn(5,10,w);`):

elif nops([args])=1 and op(1,[args])=Wt1 then
print(`Wt1(L,w,n) The simple weight of the list L . Try:`):
print(`Wt1([2,5],w,10);`):

elif nops([args])=1 and op(1,[args])=Wtpi then
print(`Wtpi(pi,w): the pure weight of the permutation pi i.e. w^(n-#FixedPoints). Try:`):
print(`Wtpi([4,2,1,3],w);`):

elif nops([args])=1 and op(1,[args])=WtpiG then
print(`WtpiG(pi,w): the GENERALIZED pure weight of the permutation pi i.e. the product of w[i] over all i that are NON-FIXED points. Try:`):
print(`WtpiG([4,2,1,3],w);`):

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




###From Findrec

 
ezraFindrec:=proc()

if args=NULL then
print(`The procedures from Findrec are : Findrec, SeqFromRec, findrec`):
print(`To get help, type ezraFindrec(procedurename);`):


elif nargs=1 and args[1]=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 nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
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 nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,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`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):


fi:
 
end:


###Findrec 
#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 (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
print(`Too many solutions`):
 RETURN(FAIL):
fi:

ope:=ope[1]:

if expand(SeqFromRec(ope,n,N,[op(1..ORDER,f)],nops(f)))<>expand(f) then
print(ope ,`did not work out`):
RETURN(FAIL):
 else
RETURN(Yafe(ope,N)[2]):
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:


YafeI:=proc(ope,yemin,N) local i,ope1,coe1,L,yemin1: 
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):
yemin1:=yemin/coe1:
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),[ope1,yemin1]:
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:





 
#SeqFromRec(ope,n,N,Ini,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
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
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)]:
gu:=expand(gu):
od:

gu:

end:


 
#SeqFromRecI(ope,yemin,n,N,Ini,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)=yemin
#extends it to the first K values.
#Try:
#SeqFromRecI(N-(n+1),n^2,N,[1],n,10);
SeqFromRecI:=proc(ope,yemin,n,N,Ini,K)
local ope1,gu,L,n1,j1,yemin1:
ope1:=YafeI(ope,yemin,N)[2]:
yemin1:=ope1[2]:
ope1:=ope1[1]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

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

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

gu:

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



###Start Inhomog.
#findrecI(f,DEGREE,ORDER,n,N): guesses an inhomogenous 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);
findrecI:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a,yemin,b:
option remember:

if (1+DEGREE)*(2+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 yemin:=0:
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:


 for j from 0 to DEGREE do
  yemin:=yemin+b[j]*n^j:
  var:=var union {b[j]}:
 od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(2+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:
 eq1:=eq1-subs(n=n0,yemin):
   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:

if nops(kv)>1 then
 print(`Too many solutions`):
 RETURN(FAIL):
fi:

ope:=subs(var1,ope):
yemin:=subs(var1,yemin):

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

ope:=normal(ope/kv[1]):
yemin:=normal(yemin/kv[1]):



if SeqFromRecI(ope,yemin,n,N,[op(1..ORDER,f)] ,nops(f))<>f then
 print([ope,yemin] , `did not work out`):
 RETURN(FAIL):
fi:


[ope,yemin]:


end:
 


#FindrecI(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 FindrecI([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
FindrecI:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

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

end:

###End From Findrec

##Start from GuessMulPol.txt

Hnm:=proc(n,m) local i:add(1/i^m,i=1..n-1):end:
#Comps1(k,r): the set of vectors of non-negative integers [v1,..., vk] that adds up to r. Try:
#Comps1(3,4);
Comps1:=proc(k,r) local gu,i,mu,mu1:
option  remember:

if k=0 then
 if r=0 then
   RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:

for i from 0 to r do
 mu:=Comps1(k-1,r-i):
 gu:=gu union {seq([op(mu1),i],mu1 in mu)}:
od:
gu:
end:


#Comps(k,r): the set of vectors of non-negative integers [v1,..., vk] that adds up to <=r. Try:
#Comps(3,4);
Comps:=proc(k,r) local r1:
{seq(op(Comps1(k,r1)),r1=0..r)}:
end:

#CompsG(A,k,r): The set of vectors of non-negative integers with k componets such that the first component is between 0 and A
#and the remaining ones have sum <=r. Try:
#CompsG(5,3,1);
CompsG:=proc(A,k,r) local a,v:
{seq(seq([a,op(v)], v in Comps(k-1,r)),a=0..A)}:
end:



#GenPol(x,k,r,a): A generic polynomial of total degree r in x[1], ..., x[k] with coefficients indexed by a,
#followed by the set of coefficients. Try:
#GenPol(x,3,3,a);
GenPol:=proc(x,k,r,a) local gu,i,j:
gu:=Comps(k,r):

[add(a[i]*mul(x[j]^gu[i][j],j=1..k),i=1..nops(gu)),{seq(a[i],i=1..nops(gu))}]:
end:


#GenPolG(x,k,A,r,a): A generic polynomial of total degree r in x[2], ..., x[k] and degree A in x[1],with coefficients indexed by a,
#followed by the set of coefficients. Try:
#GenPolG(x,3,3,1,a);
GenPolG:=proc(x,k,A,r,a) local gu,i,j:
gu:=CompsG(A,k,r):
[add(a[i]*mul(x[j]^gu[i][j],j=1..k),i=1..nops(gu)),{seq(a[i],i=1..nops(gu))}]:
end:

#ApplyPolSeq(P,x,LL): Inputs a polynomial P of x[1], ..., x[k] (where k=nops(LL)) and a list LL of length k
#consisting of the values of k sequences LL[1], ..., LL[k] all of length N, outputs
#the sequence of length N whose n-th entry is P(LL[1][n], ..., LL[k][n]). Try
#ApplyPolSeq(x[1]+x[2]+2*1]*x[2],[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]);
ApplyPolSeq:=proc(P,x,LL) local k,i,n,N:

k:=nops(LL):

if not  (type(LL,list)  and nops(LL)>0 and {seq(type(LL[i],list),i=1..k)}={true} and nops({seq(nops(LL[i]),i=1..k)})=1) then
 print(LL, `must be a non-empty list of lists all of the same length `):
 RETURN(FAIL):
fi:

N:=nops(LL[1]):
[seq(subs({seq(x[i]=LL[i][n],i=1..k)},P),n=1..N)]:
end:

#GuessPol1(LL,L1,r,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and
#a positive integer r, outputs a polynomial P of total degree <=r , if it exists, such that 
#L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Otherwise it returns FAIL. Try:
# GuessPol1([[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]], ApplyPolSeq(x[1]+x[2],x,[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]),1,x);
GuessPol1:=proc(LL,L1,r,x) local i,P,eq,var,gu,a,var1,N,k,n:
k:=nops(LL):

if not  (type(LL,list)  and nops(LL)>0 and {seq(type(LL[i],list),i=1..k)}={true} and nops({seq(nops(LL[i]),i=1..k)})=1) then
 print(LL, `must be a non-empty list of lists all of the same length `):
 RETURN(FAIL):
fi:

N:=nops(LL[1]):

if not (type(L1,list) and nops(L1)=N )then
 print(L1, `must be a list of numbers of length`, N):
 RETURN(FAIL):
fi:

gu:=GenPol(x,k,r,a):

P:=gu[1]: var:=gu[2]:

if N-nops(var)<=5 then
 print(`The length of the sequences must be at least`, nops(var)+6 ):
 RETURN(FAIL):
fi:

eq:={seq(subs({seq(x[i]=LL[i][n],i=1..k)},P)-L1[n],n=1..N)}:
var1:=solve(eq,var):

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

subs(var1,P):


end:



#GuessPol(LL,L1,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and
#a positive integer r, outputs a polynomial P such that 
#L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Try:
# GuessPol([[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]], ApplyPolSeq(x[1]+x[2],x,[[seq(i,i=1..10)], [seq(Hnm(i,1),i=1..10)]]),x);
GuessPol:=proc(LL,L1,x) local r,gu,N,i,k:
k:=nops(LL):

if not  (type(LL,list)  and nops(LL)>0 and {seq(type(LL[i],list),i=1..k)}={true} and nops({seq(nops(LL[i]),i=1..k)})=1) then
 print(LL, `must be a non-empty list of lists all of the same length `):
 RETURN(FAIL):
fi:

N:=nops(LL[1]):

if not (type(L1,list) and nops(L1)=N )then
 print(L1, `must be a list of numbers of length`, N):
 RETURN(FAIL):
fi:


for r from 0 while nops(Comps(k,r))<=N-6 do
gu:=GuessPol1(LL,L1,r,x):
if gu<>FAIL then
 RETURN(gu):
fi:
od:

FAIL:
end:




#GuessPol1G(LL,L1,A,r,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and
#a positive integers A and r, outputs a polynomial P of degree<=A in x[1] and total degree <=r in {x[2], ..., x[r]}, if it exists, such that 
#L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Otherwise it returns FAIL. Try:
# GuessPol1G([[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]], ApplyPolSeq(x[1]^3+x[2],x,[[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]]),3,1,x);
GuessPol1G:=proc(LL,L1,A,r,x) local i,P,eq,var,gu,a,var1,N,k,n:
k:=nops(LL):

if not  (type(LL,list)  and nops(LL)>0 and {seq(type(LL[i],list),i=1..k)}={true} and nops({seq(nops(LL[i]),i=1..k)})=1) then
 print(LL, `must be a non-empty list of lists all of the same length `):
 RETURN(FAIL):
fi:

N:=nops(LL[1]):

if not (type(L1,list) and nops(L1)=N )then
 print(L1, `must be a list of numbers of length`, N):
 RETURN(FAIL):
fi:

gu:=GenPolG(x,k,A,r,a):

P:=gu[1]: var:=gu[2]:

if N-nops(var)<=5 then
 print(`The length of the sequences must be at least`, nops(var)+6 ):
 RETURN(FAIL):
fi:

eq:={seq(subs({seq(x[i]=LL[i][n],i=1..k)},P)-L1[n],n=1..N)}:
var1:=solve(eq,var):

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

subs(var1,P):


end:


##End from GuessMulPol.txt



###Start From Stat


 #Moms(f,t,K): inputs a prob. gen. function and outputs the first K moments
 Moms:=proc(f,t,K) local L,i,F:

 if subs(t=1,f)=0 then
 RETURN(FAIL):
 fi:
 
 F:=f/subs(t=1,f):
 
 F:=f:

  L:=[]:

   for i from 1 to K do
    F:=t*diff(F,t):
    L:=[op(L),subs(t=1,F)]:
   od:
   L:

 end:

#MomsAM(f,t,K): inputs a prob. gen. function and outputs the first K moments about the mean. The first output is
#the expectation, the second the variance, etc.
 MomsAM:=proc(f,t,K) local L,i,F,F1,ave:
 
if subs(t=1,f)=0 then
  RETURN(FAIL):
fi:

 F:=f/subs(t=1,f):

ave:=subs(t=1,diff(F,t)):

L:=Moms(expand(F/t^ave),t,K):

[ave,op(2..K,L)]:

 end:


#Alpha(f,t,K): inputs a prob. gen. function and outputs the expectation, variance, and the scaled moments up to the K-th
#Try:
#Alpha((1+x)^100,x,6);
Alpha:=proc(f,t,K) local gu,i:

gu:= MomsAM(f,t,K):

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

evalf([op(1..2,gu),seq(gu[i]/gu[2]^(i/2),i=3..K)]):
end:

###End From Stat




#Pn(n): The set of lists [1,i1, .., ir] such that 1<=i1<i2<...<ir<=n Try: Pn(5);
Pn:=proc(n) local gu,gu1:option remember:if n=1 then  RETURN({[]}): fi: gu:=Pn(n-1): gu union {seq([op(gu1),n],gu1 in gu)}: end:


###start general case
#Wt1G(L,w,n) The general weight, in w[1], ..., w[n], of the list L
Wt1G:=proc(L,w,n) local i1: if L=[] then RETURN(1/n): else w[1]/n*mul(w[L[i1]]/(n-L[i1]+1),i1=1..nops(L)):fi: end:


#Anw1G(n,w): The full probability generating function, in w[1], ..., w[n], for n passengers, the first being absent-minded.from the definition, by brute force
#summing the weights of all possible permutations
Anw1G:=proc(n,w) local gu,gu1: gu:=Pn(n):add(Wt1G(gu1,w,n), gu1 in gu): end:

#WrnG(r,n,w): the weight enumerator, in w[1], ..., w[n] of all sequences L that end in r. Using dynamical programing
WrnG:=proc(r,n,w)  local s,gu: option remember:
if not type(r,integer) and r>=2 and r<=n then  print(`bad input`):  RETURN(FAIL): fi:

gu:=1/n/(n-r+1)*w[1]*w[r]:
for s from 2 to r-1 do  gu:=expand(gu+w[r]*WrnG(s,n,w)/(n-r+1)): od: gu:
end:


#Anw1FaG(n,w): The full probability generating function for n passengers, the first being absent-minded, using dynamical programming. Try:
#Anw1FaG(4,w);
Anw1FaG:=proc(n,w) local r: option remember: 1/n+add(WrnG(r,n,w),r=2..n): end:

#Anw1FoG(n,w): The full probability generating function for n passengers, the first being absent-minded, using the closed-form formula. Try:
#Anw1FoG(4,w);
Anw1FoG:=proc(n,w) local r: (1-w[1])/n+ w[1]/n!*mul(w[r]+n+1-r,r=2..n):end:

#Anw2FoG(n,w): The full probability generating function for n passengers, the first and second being absent-minded, using the closed-form formula. Try:
#Anw2FoG(4,w);
Anw2FoG:=proc(n,w) local r: 

if n<2 then 
 print(n, `must be at least 2`):
 RETURN(FAIL):
fi:

2/n!*w[1]*w[2]*mul(2*w[r]+n+1-r,r=3..n)+w[1]*(1-w[2])*mul(w[r]+n+1-r,r=3..n)/n!+w[2]*(1-w[1])*mul(w[r]+n+1-r,r=3..n)/n!+(1-w[1])*(1-w[2])/(n*(n-1)):

end:

###end general case

###start simple case
#Wt1(L,w,n) The simple weight of the list L . Try:
#Wt1([2,5],w,10);
Wt1:=proc(L,w,n) local i1: if L=[] then RETURN(1/n): else w^(nops(L)+1)/n*mul(1/(n-L[i1]+1),i1=1..nops(L)):fi: end:


#Anw1(n,w): The full probability generating function for n passengers, the first being absent-minded.from the definition, by brute force
#summing the weights of all possible permutations
Anw1:=proc(n,w) local gu,gu1: gu:=Pn(n):add(Wt1(gu1,w,n), gu1 in gu): end:

#Wrn(r,n,w): the weight enumerator of all sequences L that end in r. Using dynamical programing. Try:
#Wrn(5,10,w);
Wrn:=proc(r,n,w)  local s,gu: option remember:
if not type(r,integer) and r>=2 and r<=n then  print(`bad input`):  RETURN(FAIL): fi:

gu:=1/n/(n-r+1)*w^2:
for s from 2 to r-1 do  gu:=expand(gu+w*Wrn(s,n,w)/(n-r+1)): od: gu:
end:


#Anw1Fa(n,w): The full probability generating function in w, for n passengers, the first being absent-minded, using dynamical programming. Try:
#Anw1Fa(4,w);
Anw1Fa:=proc(n,w) local r: option remember: 1/n+add(Wrn(r,n,w),r=2..n): end:

#Anw1Fo(n,w): The full probability generating function in w, for n passengers, the first being absent-minded, using the closed-form formula. Try:
#Anw1Fo(4,w);
Anw1Fo:=proc(n,w) local r: (1-w)/n+ w/n!*mul(w+n+1-r,r=2..n):end:
###end simple case

#Simk(k,n): Outputs a random permutation obtained from the absent-minded passenger problem with n passengers, the first k of whom
#are absent-minded. Try:
#Simk(3,20);
Simk:=proc(k,n) local T,X,Empties,p,i:

for i from 1 to n do
 T[i]:=X:
od:

Empties:={seq(i,i=1..n)}:

for i from 1 to k do
  p:=Empties[rand(1..nops(Empties))()]:
  Empties:=Empties minus {p}:
   T[p]:=i:
 od:


for i  from k+1 to n do
 if member(i, Empties) then
  T[i]:=i:
 Empties:=Empties minus {i}:
  else
  p:=Empties[rand(1..nops(Empties))()]:
   T[p]:=i:
    Empties:=Empties minus {p}:
 fi:
od:
[seq(T[i],i=1..n)]:
end:


#SimkV(k,n): Outputs a random permutation obtained from the absent-minded passenger problem with n passengers, the first k of whom
#are absent-minded. Verbose version, prints intermediage steps. Try:
#SimkV(3,20);
SimkV:=proc(k,n) local T,X,Empties,p,i,i1:

for i from 1 to n do
 T[i]:=X:
od:

Empties:={seq(i,i=1..n)}:


print([seq(T[i1],i1=1..n)]):

for i from 1 to k do
  p:=Empties[rand(1..nops(Empties))()]:
  Empties:=Empties minus {p}:
   T[p]:=i:
 od:

print([seq(T[i],i=1..n)]):

for i  from k+1 to n do
 if member(i, Empties) then
  T[i]:=i:
 Empties:=Empties minus {i}:
  else
  p:=Empties[rand(1..nops(Empties))()]:
   T[p]:=i:
    Empties:=Empties minus {p}:
 fi:
print([seq(T[i1],i1=1..n)]):
od:
[seq(T[i1],i1=1..n)]:
end:
 



##begin simulation

#Wtpi(pi,w): the pure weight of the permutation pi i.e. w^(n-#FixedPoints). Try:
#Wtpi([4,2,1,3],w);
Wtpi:=proc(pi,w) local i,gu:
gu:=1:

for i from 1 to nops(pi) do
 if pi[i]<>i then
  gu:=gu*w:
 fi:
od:
gu:
end:
  

#AnwkSim(n,w,k,N): simulates the absent-minded passenger with k absent minded passengers and a total of n passengers 
#where the first k ones are absent-minded
#N times and outputs
#the fraction of those where pi[n]<>n and the generating function according to the number of fixed points. Try:
#AnwkSim(1,20,1000,w);
AnwkSim:=proc(n,w,k,N) local i,co,pi,gu:

co:=0:
gu:=0:

for i from 1 to N do
 pi:=Simk(k,n):
 if pi[n]=n then
  co:=co+1:
 fi:
 gu:=gu+Wtpi(pi,w):
od:

[evalf(co/N),evalf(sort(gu/N))]:
end:


#WtpiG(pi,w): The detailed pure weight of a permutation pi
#WtpiG([4,2,1,3],w);
WtpiG:=proc(pi,w) local i,gu:
gu:=1:

for i from 1 to nops(pi) do
 if pi[i]<>i then
  gu:=gu*w[i]:
 fi:
od:
gu:
end:
  


#AnwkGsim(n,w,k,N): simulates the absent-minded passenger with k absent minded passengers and a total of n passengers 
#detailed weight
#where the first k ones are absent-minded
#N times and outputs
#the fraction of those where pi[n]<>n and the generating function according to the number of non-fixed points
AnwkGsim:=proc(n,w,k,N) local i,co,pi,gu:

co:=0:
gu:=0:

for i from 1 to N do
 pi:=Simk(k,n):
 if pi[n]=n then
  co:=co+1:
 fi:
 gu:=gu+WtpiG(pi,w):
od:

[evalf(co/N),evalf(sort(gu/N))]:
end:


##end simulation

#FC1(pi,i): the cycle belonging to i in the permutation pi. Try
#FC1([2,3,1],1);
FC1:=proc(pi,i) local c,gu:
if pi[i]=i then
 RETURN([i]):
fi:

gu:=[i]:
c:=pi[i]:

while c<>i do
gu:=[op(gu),c]:
c:=pi[c]:
od:
gu:
end:

#Cyc(pi): the list of cycles of the permutation pi, except for fixed points, try:
#Cyc([3,1,2,4]);
Cyc:=proc(pi) local n,gu,mu,i,lu:
n:=nops(pi):
mu:={seq(i,i=1..n)}:

gu:=[]:

while mu<>{} do
lu:=FC1(pi,min(op(mu))):
mu:=mu minus {op(lu)}:

if nops(lu)>1 then
 gu:=[op(gu),lu]:
fi:
od:
gu:
end:


#Childn(St): Given a state St=[Pas,Seats,AM] where Pa is the current list of passengers (in order), Seats is the current list of seats (in order)
#and AM is the subset of  Pas that is abent-minded, outputs the set of children, i.e. the set of possible states
#that results from seating the first passender in Pa. Try:
#Childn([{1,2,3,4},{1,2,3,4},{1}]);
Childn:=proc(St) local Pas,Seats,AM,p,s:
Pas:=St[1]: Seats:=St[2]: AM:=St[3]:
p:=min(op(Pas)):


if not member(p,AM) and member(p,Seats) then
Pas:=Pas minus {p}:
Seats:=Seats minus {p}:
RETURN({[Pas,Seats,AM]}):
else
RETURN({seq([Pas minus {p},Seats minus {s},AM minus {p}],s in Seats)}):
fi:

end:

#ChildnWWg(St,w): Given a state St=[Pas,Seats,AM] where Pa is the current list of passengers (in order), Seats is the current list of seats (in order)
#and AM is the subset of  Pas that is abent-minded, outputs the set of 
#pairs [child, weight] i.e. the set of possible states
#that results from seating the first passender in Pa. Try:
#ChildnWWg([{1,2,3,4},{1,2,3,4},{1}],w);
ChildnWWg:=proc(St,w) local Pas,Seats,AM,p,s,gu:
Pas:=St[1]: Seats:=St[2]: AM:=St[3]:
p:=min(op(Pas)):


if not member(p,AM) and member(p,Seats) then
Pas:=Pas minus {p}:
Seats:=Seats minus {p}:
RETURN({[[Pas,Seats,AM],1]}):
else
gu:={}:
for s in Seats do
if p=s then
 gu:=gu union {[[Pas minus {p},Seats minus {s},AM minus {p}],1/nops(Seats)]}:
else
 gu:=gu union {[[Pas minus {p},Seats minus {s},AM minus {p}],w[p]/nops(Seats)]}:
fi:
od:
RETURN(gu):
fi:
end:

#ChildnWW(St,w): Given a state St=[Pas,Seats,AM] where Pa is the current list of passengers (in order), Seats is the current list of seats (in order)
#and AM is the subset of  Pas that is abent-minded, outputs the set of 
#pairs [child, weight] i.e. the set of possible states
#that results from seating the first passender in Pa. Try:
#ChildnWW([{1,2,3,4},{1,2,3,4},{1}],w);
ChildnWW:=proc(St,w) local Pas,Seats,AM,p,s,gu:
Pas:=St[1]: Seats:=St[2]: AM:=St[3]:
p:=min(op(Pas)):


if not member(p,AM) and member(p,Seats) then
Pas:=Pas minus {p}:
Seats:=Seats minus {p}:
RETURN({[[Pas,Seats,AM],1]}):
else
gu:={}:
for s in Seats do
if p=s then
 gu:=gu union {[[Pas minus {p},Seats minus {s},AM minus {p}],1/nops(Seats)]}:
else
 gu:=gu union {[[Pas minus {p},Seats minus {s},AM minus {p}],w/nops(Seats)]}:
fi:
od:
RETURN(gu):
fi:
end:

#AllStates(n,k): The set of all states coming from n passengers {1, ....,n} the first k of them being absent-minded. Try:
#AllStates(5,1);
AllStates:=proc(n,k) local cur,cur1,gu,gu1,i1:

cur:={[{seq(i1,i1=1..n)},{seq(i1,i1=1..n)},{seq(i1,i1=1..k)}]}:
gu:=cur:

cur1:= {seq(op(Childn(gu1)),gu1 in cur)}:

while cur1<>{} do
 gu:=gu union cur1:
 cur:=cur1:
 cur1:= {seq(op(Childn(gu1)),gu1 in cur)}:
od:

gu:

end:





#Fstate(St,w): the weight of  a state St. Try:
#Fstate([{1,2,3,4},{1,2,3,4},{}],w);
Fstate:=proc(St,w) local gu,i:
option remember:
if St=[{},{},{}] then
 RETURN(1):
fi:

gu:=ChildnWW(St,w):
expand(add(gu[i][2]*Fstate(gu[i][1],w),i=1..nops(gu))):

end:


#Anwk(n,w,k): the generating function in w whose coefficient of w^i is the probability that exactly i passengers will seat in 
#a different seat that has been assigned to them of there are n passengers {1, ..., n} the first k of whom are absent-minded.
#Try:Anwk(5,w,1);
Anwk:=proc(n,w,k)  local i: option rememeber: Fstate([{seq(i,i=1..n)},{seq(i,i=1..n)},{seq(i,i=1..k)}],w): end:


#FstateG(St,w): the generalized weight of  a state St. Try:
#FstateG([{1,2,3,4},{1,2,3,4},{1}],w);
FstateG:=proc(St,w) local gu,i:
option remember:
if St=[{},{},{}] then
 RETURN(1):
fi:

gu:=ChildnWWg(St,w):
expand(add(gu[i][2]*FstateG(gu[i][1],w),i=1..nops(gu))):

end:


#AnwkG(n,w,k): the generalized generating function in w[1], ..., w[n] 
#a different seat that has been assigned to them of there are n passengers {1, ..., n} the first k of whom are absent-minded.
#Try:AnwkG(5,w,1);
AnwkG:=proc(n,w,k)  local i: FstateG([{seq(i,i=1..n)},{seq(i,i=1..n)},{seq(i,i=1..k)}],w): end:


#PrS(n,k,S): Inputs positive integers n and k, and a subset S of {1,...n}, outputs the probability
#that with n passengers, the first k of whom are absent-minded, the members of S will all sit in the wrong seat. Try:
#PrS(10,1,{10});
PrS:=proc(n,k,S) local gu,S1,i,w:
gu:=AnwkG(n,w,k):
S1:={seq(i,i=1..n)} minus S:
gu:=subs({seq(w[i]=1, i in S1)},gu):

for i in S do
 gu:=coeff(gu,w[i],1):
od:

gu:

end:

#IniState(n,k): the initial state with n passengers and k absent-minded ones. Try:
#IniState(6,2);
IniState:=proc(n,k) local i:
[{seq(i,i=1..n)},{seq(i,i=1..n)},{seq(i,i=1..k)}]:
end:



#EVLL(n,k): the expectation and variance according to Norbert Henze and Guenter Last
EVLL:=proc(n,k) local i:
[
k*(1+Hnm(n,1)-Hnm(k+1,1)),
k*( (2*(1-n+k*n)-n^2-k)/n^2/(n-1)+2/(n*k)+(1+2/n)*(Hnm(n+1,1)-Hnm(k+1,1))-k*add(1/i^2,i=k+1..n))
]:
end:


#Vark(k,n,Hn): The explicit expression for the variance according to Hence and last for W_{n,k}. Try: 
#Var(2,n,Hn)
Vark:=proc(k,n,Hn):
normal(k*( (2*(1-n+k*n)-n^2-k)/n^2/(n-1)+2/(n*k)+(1+2/n)*(Hn[1]+1/n-Hnm(k,1)-1/k)-k*(Hn[2]+1/n^2-Hnm(k,2)-1/k^2))):
end:


#CoeffS(n,k,S): Inputs positive integers n and k, and a subset S={i1, ..., ik}, of {1,...n}, outputs the coefficient of
#w[i1]...w[ik] in AnwkG(n,w,k); Try:
#CoeffS(10,2,{3,10});
CoeffS:=proc(n,k,S) local gu,i,w:
gu:=AnwkG(n,w,k):

for i from 1 to n do
 if member(i,S) then
  gu:=coeff(gu,w[i],1):
 else
  gu:=coeff(gu,w[i],0):
 fi:
od:
gu:

end:








#AnwkSeq(N,w,k): The first N terms of the sequence Anwk(n,w,k). Try:
#AnwkSeq(40,w,2);
AnwkSeq:=proc(N,w,k) local n:
[seq(Anwk(n,w,k),n=1..N)]:
end:




#Opers(k,w,n,N): The linear recurrence equation annihilating the sequence Anwk(n,w,k) for k=1,2,3, starting at n=k. Try:
#Opers(3,w,n,N)
Opers:=proc(k,w,n,N) :
if k=1 then

n*(n+w)/(n+2)/(n+1)-(2*n+w+1)/(n+2)*N+N^2:

elif k=2 then

-n*(n+2*w)*(n+w)/(n+4)/(n+3)/(n+2)+(3*n^2+6*n*w+2*w^2+3*n+3*w+1)/(n+3)/(n+4)*N-3*(n+w+1)/(n+4)*N^2+N^3:

elif k=3 then
n*(n+2*w)*(n+3*w)*(n+w)/(n+5)/(n+4)/(n+3)/(n+6)-(3*w+1+2*n)*(2*n^2+6*n*w+2*w^2+2*n+3*w+1)/(n+5)/(n+4)/(n+6)*N+(6*n^2+18*n*w+11*w^2+12*n+18*w+7)/(n+6)/(n+5)*N^2-2*(2
*n+3*w+3)/(n+6)*N^3+N^4:

elif k=4 then
-n*(n+4*w)*(n+2*w)*(n+3*w)*(n+w)/(n+5)/(n+4)/(n+8)/(n+7)/(n+6)+(5*n^4+40*n^3*w+105*n^2*w^2+100*n*w^3+24*w^4+10*n^3+60*n^2*w+105*n*w^2+50*w^3+10*n^2+40*n*w+35*w^2+5*
n+10*w+1)/(n+5)/(n+8)/(n+7)/(n+6)*N-5*(2*w+1+n)*(2*n^2+8*n*w+5*w^2+4*n+8*w+3)/(n+8)/(n+7)/(n+6)*N^2+5*(2*n^2+8*n*w+7*w^2+6*n+12*w+5)/(n+7)/(n+8)*N^3-5*(n+2*w+2)/(n+
8)*N^4+N^5

else
FAIL:

fi:

end:





#AnwkSeqF(N1,w,k): The first N1 terms of the sequence Anwk(n,w,k). Fast way.Try:
#AnwkSeqF(40,w,2);
AnwkSeqF:=proc(N1,w,k) local n,ope,N,Ini,gu,i:

if not (1<=k and k<=4) then
 RETURN([seq(Anwk(n,w,k),n=1..N1)]):
fi:

ope:=Opers(k,w,n,N):
Ini:=[seq(Anwk(i,w,k),i=k..degree(ope,N)+k-1)]:
gu:=SeqFromRec(ope,n,N,Ini,N1-(k-1)):

[0$(k-1),op(gu)]:

end:


##############






####new GuessPol1G



#Comps1New(k,r): the set of vectors of non-negative integers [v1,..., vk] that adds up to r. Try:
#Comps1New(3,4);
Comps1New:=proc(k,r) local gu,i,mu,mu1:
option remember:

if r=0 then
 RETURN({[0$k]}):
fi:

if r<0 then
 RETURN({}):
fi:

if k=0 then
 if r=0 then
  RETURN({[]}):
 else
 RETURN({}):
fi:
fi:

gu:={}:

for i from 0 to trunc(r/k) do
 mu:=Comps1New(k-1,r-i*k):
 gu:=gu union {seq([op(mu1),i],mu1 in mu)}:
od:
gu:

end:


#CompsNew(k,r): the set of vectors of non-negative integers [v1,..., vk] with weight  <=r. Try:
#CompNews(3,4);
CompsNew:=proc(k,r) local r1:
{seq(op(Comps1New(k,r1)),r1=0..r)}:
end:


#CompsGnew(A,k,r): The set of vectors of non-negative integers with k componets such that the first component is between 0 and A
#and the remaining ones have weight <=r. Try:
#CompsGnew(5,3,1);
CompsGnew:=proc(A,k,r) local a,v:
{seq(seq([a,op(v)], v in CompsNew(k-1,r)),a=0..A)}:
end:

#GenPolGnew(x,k,A,r,a): A generic polynomial of total degree r in x[2], ..., x[k] and degree A in x[1],with coefficients indexed by a,
#followed by the set of coefficients. Try:
#GenPolGnew(x,3,3,1,a);
GenPolGnew:=proc(x,k,A,r,a) local gu,i,j:
gu:=CompsGnew(A,k,r):
[add(a[i]*mul(x[j]^gu[i][j],j=1..k),i=1..nops(gu)),{seq(a[i],i=1..nops(gu))}]:
end:

#GuessPol1Gnew(LL,L1,A,r,x): inputs a list of lists of the same length LL and a target list L1 of the same lengths as the lists of LL and
#a positive integers A and r, outputs a polynomial P of degree<=A in x[1] and total weight <=r in {x[2], ..., x[r]}, if it exists, such that 
#L1[n]=P(LL[1][n], ..., LL[k][n]) for n from 1 to N. Otherwise it returns FAIL. Try:
# GuessPol1Gnew([[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]], ApplyPolSeq(x[1]^3+x[2],x,[[seq(i,i=1..20)], [seq(Hnm(i,1),i=1..20)]]),3,1,x);
GuessPol1Gnew:=proc(LL,L1,A,r,x) local i,P,eq,var,gu,a,var1,N,k,n:
k:=nops(LL):

if not  (type(LL,list)  and nops(LL)>0 and {seq(type(LL[i],list),i=1..k)}={true} and nops({seq(nops(LL[i]),i=1..k)})=1) then
 print(LL, `must be a non-empty list of lists all of the same length `):
 RETURN(FAIL):
fi:

N:=nops(LL[1]):

if not (type(L1,list) and nops(L1)=N )then
 print(L1, `must be a list of numbers of length`, N):
 RETURN(FAIL):
fi:

gu:=GenPolGnew(x,k,A,r,a):

P:=gu[1]: var:=gu[2]:

if N-nops(var)<=5 then
 print(`The length of the sequences must be at least`, nops(var)+6 ):
 RETURN(FAIL):
fi:

eq:={seq(subs({seq(x[i]=LL[i][n],i=1..k)},P)-L1[n],n=1..N)}:
var1:=solve(eq,var):

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

subs(var1,P):


end:






#MomPask1(n,Hn,r): Guesses an expression for the r-th moment about the  mean of W_{n,1}. Try:
#MomPask1(n,Hn,4);
MomPask1:=proc(n,Hn,r) local N,vu,w,r1,i,x,gu,lu:

N:=nops(CompsGnew(1,r,r))+10:

vu:=AnwkSeqF(1+N,w,1):
gu:=[seq(i*MomsAM(vu[i],w,r)[r],i=1..1+N)]:

lu:=GuessPol1Gnew([[seq(i,i=1..1+N)],seq([seq(Hnm(i,r1),i=1..1+N)],r1=1..r)],gu,1,r,x):
normal(subs({x[1]=n,seq(x[r1+1]=Hn[r1],r1=1..r)},lu)/n):
end:


#MomPas(k,n,Hn,r): Finds an explicit expression for the expected number (if r=1)
#and for the r-th moment about the mean for r>=2 of displaced passengers with k absent-minded passengers. Try:,
#in terms of n and the Harmonic numbers Hn. Try:
#MomPas(1,n,Hn,1);
MomPas:=proc(k,n,Hn,r) local gu,i,w,lu,x,N,r1,vu:

if not (type(r,integer) and r>=1 and r<=4) then
 print(r, `must be an integer between 1 and 4`):
 RETURN(FAIL):
fi:

if r=1 then
N:=10:

vu:=AnwkSeqF(k+N,w,k):

gu:=[seq(MomsAM(vu[i],w,r)[r],i=k..k+N)]:
lu:=GuessPol1Gnew([[seq(i,i=k..k+N)],seq([seq(Hnm(i,r1),i=k..k+N)],r1=1..r)],gu,0,1,x):

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

RETURN(expand(subs({x[1]=n,seq(x[r1+1]=Hn[r1],r1=1..r)},lu))):
fi:


if r=2 then
N:=18:
vu:=AnwkSeqF(k+N,w,k):
gu:=[seq(i*(i-1)*MomsAM(vu[i],w,r)[r],i=k..k+N)]:
lu:=GuessPol1Gnew([[seq(i,i=k..k+N)],seq([seq(Hnm(i,r1),i=k..k+N)],r1=1..r)],gu,2,2,x):

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

RETURN(normal(subs({x[1]=n,seq(x[r1+1]=Hn[r1],r1=1..r)},lu)/(n*(n-1)))):
fi:



if r=3 then


if k=1 then

N:=20:


vu:=AnwkSeqF(k+N,w,k):

gu:=[seq(i*MomsAM(vu[i],w,r)[r],i=k..k+N)]:

lu:=GuessPol1Gnew([[seq(i,i=k..k+N)],seq([seq(Hnm(i,r1),i=k..k+N)],r1=1..r)],gu,1,3,x):
if lu=FAIL then
  RETURN(FAIL):
fi:


RETURN(normal(subs({x[1]=n,seq(x[r1+1]=Hn[r1],r1=1..r)},lu)/n)):

elif k=2 then
N:=26:

vu:=AnwkSeqF(k+N,w,k):
gu:=[seq(i*(i-1)*MomsAM(vu[i],w,r)[r],i=k..k+N)]:
lu:=GuessPol1Gnew([[seq(i,i=k..k+N)],seq([seq(Hnm(i,r1),i=k..k+N)],r1=1..r)],gu,2,3,x):



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

RETURN(normal(subs({x[1]=n,seq(x[r1+1]=Hn[r1],r1=1..r)},lu)/n/(n-1) )):

else
N:=100:

vu:=AnwkSeqF(k+N,w,k):

gu:=[seq(i*(i-1)*MomsAM(vu[i],w,r)[r],i=k..k+N)]:

lu:=GuessPol1Gnew([[seq(i,i=k..k+N)],seq([seq(Hnm(i,r1),i=k..k+N)],r1=1..r)],gu,2,3,x):

RETURN(normal(subs({x[1]=n,seq(x[r1+1]=Hn[r1],r1=1..r)},lu)/(n*(n-1)))):


fi:

fi:

if r=4 then


if k=1 then
N:=30:
vu:=AnwkSeqF(k+N,w,k):
gu:=[seq(i*MomsAM(vu[i],w,r)[r],i=k..k+N)]:

lu:=GuessPol1Gnew([[seq(i,i=k..k+N)],seq([seq(Hnm(i,r1),i=k..k+N)],r1=1..r)],gu,1,4,x):

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

RETURN(normal(subs({x[1]=n,seq(x[r1+1]=Hn[r1],r1=1..r)},lu)/n)):


else 
N:=50:

vu:=AnwkSeqF(k+N,w,k):
gu:=[seq(i*(i-1)*MomsAM(vu[i],w,r)[r],i=k..k+N)]:
lu:=GuessPol1Gnew([[seq(i,i=k..k+N)],seq([seq(Hnm(i,r1),i=k..k+N)],r1=1..r)],gu,2,4,x):



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


RETURN(normal(subs({x[1]=n,seq(x[r1+1]=Hn[r1],r1=1..r)},lu)/n/(n-1) )):

fi:


fi:



end:





#MomPask1Asy(n,r,K): The asymptotic expression for  the r-th moment about the  mean of W_{n,1} up to the K-th order. Try:
#MomPask1Asy(n,4,6);
MomPask1Asy:=proc(n,r,K) local gu,r1,i,Hn:
gu:=MomPask1(n,Hn,r):

gu:=subs({seq(Hn[r1]=sum(1/i^r1,i=1..n-1),r1=1..r)},gu):
asympt(gu,n,K):

end:


#AlphaPask1Asy(n,r,K): The asymptotic expression for  the scaled r-th moment about the  mean of W_{n,1} up to K terms
#if r is even, and its square it r is odd.
#AlphaPask1Asy(n,4,6);
AlphaPask1Asy:=proc(n,r,K) local gu2,gu, i,Hn,r1:

gu:=MomPask1(n,Hn,r):
gu2:=MomPask1(n,Hn,2):

if r mod 2=0 then
 gu:=gu/gu2^(r/2):
else
 gu:=gu^2/gu2^r:
fi:
gu:=subs({seq(Hn[r1]=sum(1/i^r1,i=1..n-1),r1=1..r)},gu):
asympt(gu,n,K):

end:


#Paper1(w,n): inputs symbols w and n and outputs an article with recurrences for
# the probability generating function for F_k(w,n+k-1): the probability generating function
#according to the number of passengers sitting in a wrong seat if there are n+k-1 passengers, k of whom are absent-minded, for
#k from 1 to 4. Try:
#Paper1(w,n):
Paper1:=proc(w,n) local k,ope,N,A,Ini,i:

print(`Linear Recurrences with Polynomial Coefficients satisfied by the Probability Generating Function According to the number of`):
print(`Passengers sitting in the wrong seat where there there are n passengers with the FIRST one being absent-minded.`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
for k from 1 to 4 do
print(``):
print(`Theorem Number`, k):

 ope:=Opers(k,w,n,N):
print(` Let `, A[k](n)(w), ` be the the  probability generating function according to the random variable NumberOfPassengersSittingInTheWrong Seat`):
print(`where there are`, n+k-1, `passengers,the first`, k, `of whom are absent-minded `):
print(`Then we have the following recurrence `):
print(``):
print(add(coeff(ope,N,i)*A[k](n+i),i=0..degree(ope,N))=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(add(coeff(ope,N,i)*A[k](n+i),i=0..degree(ope,N))=0):
print(``):
Ini:=[seq(Anwk(i,w,k),i=k..degree(ope,N)+k-1)]:
print(``):
print(`subject to the initial conditions`):
print(``):
print(seq(A[k](i)(w)=Ini[i],i=1..nops(Ini))):
print(``):
print(``):
print(`and in Maple notation`):
print(``):
lprint(seq(A[k](i)(w)=Ini[i],i=1..nops(Ini))):
print(``):
od:

end:

#Paper2(n,Hn,M1,M2): inputs symbols n, H, and positive integers M1, and M2, outputs explicit expressions for
#the r-th moment about the mean, for r between 1 and M1, of the random variable "Number of Passengers seating in the wrong seat" where
#the FIRST passenger is absent-minded, as well as asymptotic expressions up to the M2-th order. Here
#Hn[r] is short for sum(1/i^r,i=1..n-1). Try:
#Paper2(n,Hn,4,6);
Paper2:=proc(n,Hn,M1,M2) local r,r1,i,gu,lu:
print(`Explicit Expressions for the Expectation, Variance and the Moments about the mean up to the`, M1, ` -th `):
print(`As well as their asymptotic to order`, M2):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

print(`Let `):
print(``):
print(Hn[r](n)=Sum(1/i^r,i=1..n-1)):
print(``):
print(`In particular, Hn[1](n) is the (n-1)-th Harmonic number.`):
print(``):
print(`Note that Hn[r](n) is the partial sum, up to n-1, of Zeta[r] `):
print(``):
print(`Consider the absent-minded passenger problem, where there are n passengers, the FIRST of whom is absent-minded `):
print(`Let X be the random variable: Number of Passengers sitting in the wrong seat`):
print(`We have the following`, M1, `theorems each of them with its own corollary `):

print(`Theorem 1: The Expectation of X is`):
print(``):
gu:=MomPask1(n,Hn,1):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):

print(``):
print(`Corollary 1: The asymptotics of the expectation to order`, M2, ` is `):
print(``):

lu:=MomPask1Asy(n,1,M2):

print(lu):
print(``):


print(`Theorem 2: The Variance of X is`):
print(``):
gu:=MomPask1(n,Hn,2):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
print(`Corollary 2: The asymptotics of the variance to order`, M2, ` is `):
print(``):
lu:=MomPask1Asy(n,2,M2):

print(lu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu):
print(``):



for r1 from 3 to M1 do


print(`Theorem Number`, r1 , `: The `, r1 , `-th moment about the mean is `):
print(``):
gu:=MomPask1(n,Hn,r1):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
print(`Corollary Number`, r1, `: The asymptotics of the variance to order`, M2, ` is `):
print(``):
lu:=MomPask1Asy(n,r1,M2):

print(lu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu):
print(``):

od:

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

end:




#Mish(pi,w): w^NumberOfNotFixedPoints
Mish:=proc(pi,w) local i,gu:
gu:=1:
for i from 1 to nops(pi) do
 if pi[i]<>i then
  gu:=gu*w:
 fi:
od:
gu:
end:

#Paper3(k,n,N,R): compares the statistical data up to the M1-th moment the result of simulating the Absent-Minded Passenger
#problem with the first k passengers being absent-minded with the actual thing using the exact probablity generating function
#Anwk(N1,w,k). Try:
#Paper3(4,30,1000,4);
Paper3:=proc(k,n,N,R) local gu,pi,w,i,lu:
gu:=0:
for i from 1 to N do
pi:=Simk(k,n):
gu:=gu+Mish(pi,w):
od:
gu:=gu/n:

lu:=AnwkSeqF(n,w,k)[n]:

[Alpha(gu,w,R), Alpha(lu,w,R)]:

print(`Simulating the Absent-Minded Passenger Problem with`, n, `passengers, the first`, k, `of whom are absent-minded`, N, `times `):
print(`Versus the exact results `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`The Theoretical historgram is`):
print(``):
print(plot([seq([i,coeff(lu,w,i)],i=0..degree(lu,w))])):
print(``):
print(`The simulated historgram is`):
print(``):
print(plot([seq([i,coeff(gu,w,i)],i=0..degree(lu,w))])):
print(``):
print(`The true average, variance, up to the first`, R, `scaled moments are `):
print(``):
print( Alpha(lu,w,R)):
print(``):
print(`The simulated average, variance, up to the first`, R, `scaled moments are `):
print(``):
print( Alpha(gu,w,R)):
print(``):


end:



#DerG(F,w,S): inputs a polynomial F in w[1], ..., w[n] (with n>=k, does not matter, n it is not part of the input) and a subset S of {1,..k}
#outputs subs({w[i]=1, in in S},F) and
#F-mul(w[i], in in S)* subs({w[i]=1, in in S},F)
DerG:=proc(F,w,S) local gu,i,lu:
gu:=subs({seq(w[i]=1,i in S)},F):
lu:=expand(F-mul(w[i],i in S)*gu):
[factor(gu),lu]:
end:

#PeelG1(w,n,k,SeqS): Inputs a symbols w, and positive integers n and a sequence of subsets of {1,...k}, outputs what happens to
#AnwkG(n,w,k) when you keep "peeling it", by applying DerG(F,w,SeqS[i]) to the list SeqS. Try
#Try: PeelG1(w,7,2,[{1,2},{2},{1}]);
PeelG1:=proc(w,n,k,SeqS) local gu,lu,F0,F1,i,ku,ru,i1:
gu:=[]:

F0:=AnwkG(n,w,k):
F1:=F0:

for i from 1 to nops(SeqS) do
lu:=DerG(F1,w,SeqS[i]):
gu:=[op(gu),lu[1]]:
F1:=lu[2]:
od:
lu:=[op(gu), factor(F1)]:

ku:=lu[nops(lu)]:

for i from nops(SeqS)  by -1 to 1 do
 ru:=SeqS[i]:
 ku:=ku+mul(w[i1],i1 in ru)*lu[i]:
od:

if expand(ku-F0)<>0 then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

lu,ku:
end:


#PeelG(w,n,k): Finds the coefficients of mul(w[i],i in S)*mul(1-w[i],i not in S) for all subsets S of {1, ..., k}  of
#AnwkG(w,n,k).  Try:
#PeelG(w,7,2);
PeelG:=proc(w,n,k) local SeqS,i,lu,i1:
SeqS:=[]:

for i from k by -1 to 1 do
lu:=choose({seq(i1,i1=1..k)},i):
SeqS:=[op(SeqS),op(lu)]:
od:


PeelG1(w,n,k,SeqS):

end:

#erk(r,k,w): e_r(w_1,...,w_k) the weight enumerator of subsets of {1,...k} with r elements where the weight is
#mul(1-w[i],i in S)*mul(w[i],i in not S). Try:
#erk(2,4,w);
erk:=proc(r,k,w) local X,i:
coeff(mul((1-w[i])*X+w[i],i=1..k),X,r):
end:



#AnwkGE(n,w,k): AnwkG(n,w,k) using the EXPLICIT formula
AnwkGE:=proc(n,w,k)   local r,j:

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

add((k-r)!/n!*erk(r,k,w)*mul((k-r)*w[j]+n+1-j,j=k+1..n),r=0..k):
end:


#AnwkE(n,w,k): Anwk(n,w,k) using the EXPLICIT formula
AnwkE:=proc(n,w,k)   local r,j:

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

add(k!/r!/n!*(1-w)^r*w^(k-r)*mul((k-r)*w+n+1-j,j=k+1..n),r=0..k):
end: