###################################################################### ##FELLER: Save this file as FELLER To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read FELLER: # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Oct./Nov. 2006 print(`Created: Oct./Nov. 2006.`): print(` This is FELLER, a Maple package `): print(`to study coin-tossing a fair coin n times`): print(`It accompanies the paper `): print(`Fully AUTOMATED Computerized Redux of Feller's (v.1) Ch. III (and Much More!)`): print(`by Doron Zeilberger`): print(`Published in the Personal Jouranl of Ekhad and Zeilberger`): print(`http://www.math.rutgers.edu/~zeilberg/pj.html `): print(`and also available from his website`): lprint(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): lprint(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg .`): print(`For a list of the procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): with(combinat): ezra1:=proc() if args=NULL then print(` The supporting procedures are: CheckFeller, CheckChungFeller`): print(`ChungFellerBrute, CollectData, FellerBruteE`): print( ` GR1, GuessY , ShowUp`): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: ChungFeller, Feller `): print(` FellerStartEreduced `): print(` ThIII41, ThIII42 , ThIII51, ThIII9 `): elif nops([args])=1 and op(1,[args])=CheckChungFeller then print(`CheckChungFeller(N): checks the ChungFeller g.f.`): print(`by comparing it with the brute force up to the Nth term`): print(`CheckChungFeller(6);`): elif nops([args])=1 and op(1,[args])=CheckFeller then print(`CheckFeller(N): checks the Feller g.f.`): print(`by comparing it with the brute force up to the Nth term`): print(`CheckFeller(6);`): elif nops([args])=1 and op(1,[args])=ChungFeller then print(`ChungFeller(z,t1,t2,t3,t4): The Chung-Feller formal-power series`): print(`whose coeff. of z^n is the weight-enumerator of all`): print(`sequences of {1,-1} of length n with the weight`): print(`weight(Sequenec):=t1^(#losing times)*t2^(#returns to the origin)`): print(`*t3^(last visit to the origin)*t4^(#sign-changes))`): print(`Do: ChungFeller(z,t1,t2,t3,t4);`): elif nops([args])=1 and op(1,[args])=ChungFellerBrute then print(`ChungFellerBrute(z,t1,t2,t3,t4,N): The first N+1 terms of`): print(`The Chung-Feller formal-power series`): print(`whose coeff. of z^n is the weight-enumerator of all`): print(`sequences of {1,-1} of length n with the weight`): print(`weight(Sequenec):=t1^(#losing times)*t2^(#returns to the origin)`): print(`*t3^(last visit to the origin)*t4^(#sign-changes))`): print(`Do: ChungFellerBrute(z,t1,t2,t3,t4);`): elif nops([args])=1 and op(1,[args])=CollectData then print(`CollectData(VarList,L,n,K1,K2): Given a list of polynomials`): print(`in the variable list VarList collects data for`): print(`the conjectured rational function `): print(`coeff(L[n+1],VarList^a)/coeff(L[n],VarList^a) `): print(`together with the n0 that it starts to be applicable to`): print(`it collects data up to the K1 entry of L and uses everything`): print(`up to the K2-th place`): print(`For example, try:`): print(`CollectData([t1,t2,t3,t4],L,n,5,25);`): elif nops([args])=1 and op(1,[args])=Feller then print(`Feller(z,t1,t2,t3,t4): The Feller formal-power series`): print(`whose coeff. of z^n is the weight-enumerator of all`): print(`sequences of {1,-1} of length n with the weight`): print(`weight(Sequenec):=t1^(#losing times)*t2^(#returns to the origin)`): print(`*t3^(last visit to the origin)*t4^(#sign-changes))`): print(`Do: Feller(z,t1,t2,t3,t4);`): elif nops([args])=1 and op(1,[args])=FellerBruteE then print(`FellerBruteE(z,t1,t2,t3,t4,N): the exapnsion of the first`): print(`N terms of the Feller generating function by brute force`): print(`For example, try:`): print(`FellerBruteE(z,t1,t2,t3,t4,3);`): elif nops([args])=1 and op(1,[args])=FellerStartEreduced then print(`FellerStartEreduced(z,t1,t2,t3,t4,N): The expansion `): print(`of the Even part of the`): print(`Feller generating function but z^2->z t1^2->t1 t3^2->t3`): print(`For example, do:`): print(`FellerStartEreduced(z,t1,t2,t3,t4,4);`): elif nops([args])=1 and op(1,[args])=GR1 then print(`GR1(S1,x): guesses a polynomial from the data-set S1`): print(`For example try: GR1({seq([i,(i+1)/(i+2)],i=1..10)},x);`): elif nops([args])=1 and op(1,[args])=ShowUp1 then print(`ShowUp1(VarList,L,i0): the exponents that show up in`): print(`the list L starting at the place i0, and remianing for ever after`): print(`ShowUp1([t1,t2,t3,t4],Feller30,5); `): elif nops([args])=1 and op(1,[args])=ShowUp then print(`ShowUp(VarList,L,i0): the exponents that show up in`): print(`the list L starting at the place i0, for the FIRST time`): print(`and remianing for ever after`): print(`ShowUp([t1,t2,t3,t4],Feller30,5); `): elif nops([args])=1 and op(1,[args])=ThIII41 then print(`ThIII41(): Proves (rigorously!) Th. III.4.1 in Feller v.1 (3rd ed)`): print(`Do: ThIII41();`): elif nops([args])=1 and op(1,[args])=ThIII42 then print(`ThIII42(): Proves (rigorously!) Th. III.4.2 in Feller v.1 (3rd ed)`): print(`Do: ThIII42();`): elif nops([args])=1 and op(1,[args])=ThIII51 then print(`ThIII51(n0): Proves (empirically) Th. III.5.1. in Feller v.1 (3rd ed)`): print(`up to n0 coin tosses: For example do`): print(`ThIII51(20);`): elif nops([args])=1 and op(1,[args])=ThIII9 then print(`ThIII9(): Proves (rigorously!) Th. III.9. in Feller v.1 (3rd ed)`): print(`due to Chung-Feller`): print(`Do: ThIII9();`): else print(`There is no ezra for`,args): fi: end: ###From GessRat######### with(SolveTools): #Bdok1(DS1,P,x): checks that P in the single variable x indeed fits the data set DS1 Bdok1:=proc(DaS,P,x) local i: for i from 1 to nops(DaS) do if subs(x=DaS[i][1],denom(P))<>0 and normal(subs(x=DaS[i][1],P)-DaS[i][2])<>0 then RETURN(false): fi: od: true: end: #GenRat1(x,dxt,dxb,a): a generic (monic top) #rational function of degree of top dxt #and degree of bottom dxb #with coeffs.parametrized by a[i] followed by the set #of coeffs. GenRat1:=proc(x,dxt,dxb,a) local i: (add(a[i]*x^i,i=0..dxt))/ add(a[i+dxt+1]*x^i,i=0..dxb) ,{seq(a[i],i=0..dxt+dxb+1)}: end: #GR1d1d2(S1,x,d1,d2): guesses a rational function in x of degree d1/d2 #that fits the data in DS1 GR1d1d2:=proc(S1,x,d1,d2) local eq,var,P,a,var1,i: P:=GenRat1(x,d1,d2,a): var:=P[2]: P:=P[1]: if nops(S1)<=d1+d2+3 then # print(`Insufficient data`): RETURN(FAIL): fi: eq:={seq(numer(normal(subs(x=S1[i][1],P)-S1[i][2])),i=1..d1+d2+3)}: var1:=Linear(eq,var): if expand(subs(var1,numer(P)))=0 then RETURN(FAIL): fi: P:=subs(var1,P): if Bdok1(S1,P,x) then RETURN(factor(normal(P))): else RETURN(FAIL): fi: end: #GR1(S1,x): guesses a polynomial from the data-set S1 #For example try: GR1({seq([i,(i+1)/(i+2)],i=1..10)},x); GR1:=proc(S1,x) local d1,d2,gu,L,i1: if {seq(S1[i1][2],i1=1..nops(S1))}={0} then RETURN(0): fi: for L from 0 to nops(S1)-4 do for d1 from 0 to L do d2:=L-d1: gu:=GR1d1d2(S1,x,d1,d2): if gu<>FAIL then RETURN(gu): fi: od: od: FAIL: end: ###End From GessRat######### #t1: #losing times #t2: #returns to the origin #t3: last visit to the origin #t4: #sign-changes Feller2:=subs({t1=1,t3=1,t4=1},Feller); Feller2e:=normal((Feller2+subs(z=-z,Feller2))/2); ChungFeller2:=subs({t1=1,t3=1,t4=1},ChungFeller); ChungFeller2e:=normal((Feller2+subs(z=-z,ChungFeller2))/2): ChungFeller3:=subs({t1=1,t2=1,t4=1},ChungFeller); Feller3e:=normal((Feller3+subs(z=-z,ChungFeller3))/2); CheckLastVisit:=simplify(Feller3e-1/sqrt(1-4*z^2)/sqrt(1-4*z^2*t3^2)): ChungFeller4:=subs({t1=1,t2=1,t3=1},ChungFeller); ChungFeller4o:=normal((Feller4-subs(z=-z,Feller4))/2)/2: #Feller(z,t1,t2,t3,t4): The Feller formal-power series #whose coeff. of z^n is the weight-enumerator of all #sequences of {1,-1} of length n with the weight #weight(Sequenec):=t1^(#losing times)*t2^(#returns to the origin) #*t3^(last visit to the origin)*t4^(#sign-changes) #Do: Feller(z,t1,t2,t3,t4); Feller:=proc(z,t1,t2,t3,t4) local G,P1,N1,P,N,psi,alpha,beta,phi: G:=normal( 1+P1+N1+( P*(1+P1+t4*N1) +N*(1+t4*P1+N1) +t4*P*N*(1+t4*P1+N1) +t4*P*N*(1+P1+t4*N1))/(1-t4^2*P*N)): psi:=(1-sqrt(1-4*z^2))/2; phi:=(1-sqrt(1-4*z^2))/2/z^2-1; alpha:=z/(1-z-psi); beta:=z/(1-z-psi): normal(subs({ P=t2*subs(z=z*t3,psi)/(1-t2*subs(z=z*t3,psi)), N=t2*subs(z=z*t1*t3,psi)/(1-t2*subs(z=z*t1*t3,psi)), P1=alpha, N1=subs(z=t1*z,beta)},G)): end: #ChungFeller(z,t1,t2,t3,t4): The Chung-Feller formal-power series #whose coeff. of z^n is the weight-enumerator of all #sequences of {1,-1} of length n that sum to 0, with the weight #weight(Sequenec):=t1^(#losing times)*t2^(#returns to the origin) #*t3^(last visit to the origin)*t4^(#sign-changes) #Do: ChungFeller(z,t1,t2,t3,t4); ChungFeller:=proc(z,t1,t2,t3,t4) local F,P1,N1,P,N,psi: F:=normal(1+(P+N+2*t4*P*N)/(1-t4^2*P*N)): psi:=(1-sqrt(1-4*z^2))/2; normal(subs({P=t2*subs(z=z*t3,psi)/(1-t2*subs(z=z*t3,psi)), N=t2*subs(z=z*t1*t3,psi)/(1-t2*subs(z=z*t1*t3,psi)) },F)): end: #FellerV(z,t1,t2,t3,t4): #Verbose version of Feller(z,t1,t2,t3,t4): FellerV:=proc(z,t1,t2,t3,t4): print(`The Feller Generating Function, which is the weight-enumerator`): print(`of all sequenecs of {-1,1} with the weight`): print(`z^(length)*t1^(#losing-times)*t2^(#returns to the origin)*`): print(`*t3^(last visit to the origin)*t4^(#sign-changes) `): print(` is: `): print(Feller(z,t1,t2,t3,t4)): end: #ThIII42(): Proves (rigorously!) Th. III.4.2. in Feller v.1 (3rd ed) #Do: ThIII42(); ThIII42:=proc() local gu,z,t1,t2,t3,t4,gue: gu:=subs({t2=1,t3=1,t4=1},Feller(z,t1,t2,t3,t4)); gue:=normal((gu+subs(z=-z,gu))/2); print(`Theorem III.4.2 in Feller vol. 1 (due to Eric Sparre-Andersen) is`): print(evalb(0=simplify(gue-1/sqrt(1-4*z^2)/sqrt(1-4*z^2*t1^2)))); end: #ThIII41(): Proves (rigorously!) Th. III.4.1. in Feller v.1 (3rd ed) #Do: ThIII41(); ThIII41:=proc() local gu,z,t1,t2,t3,t4,gue: gu:=subs({t1=1,t2=1,t4=1},Feller(z,t1,t2,t3,t4)); gue:=normal((gu+subs(z=-z,gu))/2); print(`Theorem III.4.1 in Feller vol. 1 is `): print(evalb(simplify(gue-1/sqrt(1-4*z^2)/sqrt(1-4*z^2*t3^2))=0)); end: #ThIII51(n0): Proves (empirically) Th. III.5.1. in Feller v.1 (3rd ed) #up to n0 coin tosses: For example do #ThIII51(20); ThIII51:=proc(n0) local n,r,Feller4,Feller4o,t1,t2,t3,t4,z: Feller4:=subs({t1=1,t2=1,t3=1},Feller(z,t1,t2,t3,t4)); Feller4o:=normal((Feller4-subs(z=-z,Feller4))/2)/2; print(`Theorem III.5.1 in Feller vol. 1 for n up to `, n0 , `is `): evalb({seq( seq (coeff(coeff(expand(taylor(Feller4o,z=0,2*n0+6)),z,2*n+1),t4,r) -binomial(2*n+1,n+r+1), r=0..n),n=0..n0)}={0}); end: #ThIII9(): Proves (rigorously!) Th. III.9. in Feller v.1 (3rd ed) #due to Chung-Feller #Do: ThIII9(); ThIII9:=proc() local gu,z,t1,t2,t3,t4,phi: phi:=(1-sqrt(1-4*z^2))/2/z^2-1; gu:=subs({t2=1,t3=1,t4=1},ChungFeller(z,t1,t2,t3,t4)); print(`Theorem III.9 in Feller vol. 1 (due to Feller-Chung, PNAS (1949)) is`): print(evalb(normal(gu-((1+phi)-t1^2*subs(z=t1*z,1+phi)) /(1-t1^2))=0)): end: #FellerStartE(z,t1,t2,t3,t4,N): The expansion of the Even part of the #Feller generating function up to z^(2*N) #For example, do: #FellerStartE(z,t1,t2,t3,t4,4); FellerStartE:=proc(z,t1,t2,t3,t4,N) local gu,i: gu:=Feller(z,t1,t2,t3,t4): gu:=normal(simplify((subs(z=-z,gu)+gu)/2)): gu:=expand(taylor(gu,z=0,2*N+6)): add(coeff(gu,z,2*i)*z^(2*i),i=0..N): end: #ChungFellerStart(z,t1,t2,t3,t4,N): The expansion of the #Chung-Feller #generating function up to z^(2*N) #For example, do: #ChungFellerStart(z,t1,t2,t3,t4,4); ChungFellerStart:=proc(z,t1,t2,t3,t4,N) local gu,i: gu:=ChungFeller(z,t1,t2,t3,t4): gu:=expand(taylor(gu,z=0,2*N+6)): add(coeff(gu,z,2*i)*z^(2*i),i=0..N): end: #ShowUp1(VarList,L,n,i0): the exponents that show up in #the list L at the place i0 and stay for ever after #Do: ShowUp1([t1,t2,t3,t4],Feller30,n,5) ; ShowUp1:=proc(VarList,L,i0) local mu,i,i1,lu,gu,lu1,j: option remember: if i0<1 then RETURN({}): fi: mu:=L[i0]: lu:= {seq([seq(degree(op(i,mu),VarList[i1]),i1=1..nops(VarList))],i=1..nops(mu))}: gu:={}: for i from 1 to nops(lu) do lu1:=op(i,lu): if {seq( evalb(Mekadem(L[j],VarList,lu1)<>0) ,j=i0+1..nops(L))}={true} then gu:=gu union {lu1}: fi: od: gu: end: #ShowUp(VarList,L,n,i0): the exponents that show up in #the list L at the place i0 for the FIRST time and stay for ever after #Do: ShowUp([t1,t2,t3,t4],Feller30,n,5) ; ShowUp:=proc(VarList,L,i0) ShowUp1(VarList,L,i0) minus ShowUp1(VarList,L,i0-1): end: #CollectData(VarList,L,n,K1,K2): Given a list of polynomials #in the variable list VarList collects data for #the conjectured rational function #coeff(L[n+1],VarList^a)/coeff(L[n],VarList^a) #together with the n0 that it starts to be applicable to #it collects data up to the K1 entry of L and uses everything #up to the K2-th place #For example, try: #CollectData([t1,t2,t3,t4],L,n,5,25); CollectData:=proc(VarList,L,n,K1,K2) local lu,i0,gu,j,mu,ku: gu:={}: for i0 from 1 to K1 do ku:=[op(1..min(nops(L),K2),L)]: lu:=ShowUp(VarList,ku,i0): for j from 1 to nops(lu) do mu:=GuessY(VarList,lu[j],L,n,i0): if not type(mu,set) then gu:=gu union {[lu[j],i0,mu]}: fi: od: od: gu: end: Hal:=proc(a1,a2,a3,a4) local lu,n0: n0:=30: lu:= Mekadem(Feller30[n0+1],[t1,t2,t3,t4],[a1,a2,a3,a4])/ binomial(2*n0-2*a3-1,n0-a3); lu: end: a1:=proc(w) local gu,s,i: gu:=0: s:=0: for i from 1 to nops(w) do s:=s+w[i]: if s<0 or s=0 and w[i]=1 then gu:=gu+1: fi: od: gu: end: a2:=proc(w) local gu,s,i: gu:=0: s:=0: for i from 1 to nops(w) do s:=s+w[i]: if s=0 then gu:=gu+1: fi: od: gu: end: a3:=proc(w) local i,j: for i from nops(w) by -1 to 0 do if add(w[j],j=1..i)=0 then RETURN(i): fi: od: end: a4:=proc(w) local gu,s,i: gu:=0: s:=0: for i from 1 to nops(w)-1 do s:=s+w[i]: if s=0 and w[i]*w[i+1]>0 then gu:=gu+1: fi: od: gu: end: #Words(n): all words of length n of {1,-1} Words:=proc(n) local gu,i: option remember: if n=0 then RETURN({[]}): fi: gu:=Words(n-1): {seq([op(gu[i]),-1],i=1..nops(gu)),seq([op(gu[i]),1],i=1..nops(gu))}: end: #CFwords(n): all the words of length n in {1,-1} that add-up to 0 CFwords:=proc(n) local gu,i,mu: mu:=Words(n): gu:={}: for i from 1 to nops(mu) do if convert(mu[i],`+`)=0 then gu:=gu union {mu[i]}: fi: od: gu: end: #FellerBruteE(z,t1,t2,t3,t4,N): the exapnsion of the first #N terms of the Feller generating function by brute force #For example, try: #FellerBruteE(z,t1,t2,t3,t4,3); FellerBruteE:=proc(z,t1,t2,t3,t4,N) local gu,n,gu1,mu1,w: gu:=0: for n from 0 to N do mu1:=Words(2*n): gu1:=0: for w in mu1 do gu1:=gu1+t1^a1(w)*t2^a2(w)*t3^a3(w)*t4^a4(w): od: gu:=gu+gu1*z^(2*n): od: gu: end: #ChungFellerBrute(z,t1,t2,t3,t4,N): the exapnsion of the first #N terms of the Chung-Feller generating function #For example, try: #ChungFellerBrute(z,t1,t2,t3,t4,3); ChungFellerBrute:=proc(z,t1,t2,t3,t4,N) local gu,n,gu1,mu1,w: gu:=0: for n from 0 to N do mu1:=CFwords(2*n): gu1:=0: for w in mu1 do gu1:=gu1+t1^a1(w)*t2^a2(w)*t3^a3(w)*t4^a4(w): od: gu:=gu+gu1*z^(2*n): od: gu: end: CheckChungFeller:=proc(N) local z,t1,t2,t3,t4: evalb(ChungFellerBrute(z,t1,t2,t3,t4,N)-ChungFellerStart(z,t1,t2,t3,t4,N)=0): end: CheckFeller:=proc(N) local z,t1,t2,t3,t4: evalb(FellerBruteE(z,t1,t2,t3,t4,N)-FellerStartE(z,t1,t2,t3,t4,N)=0): end: #FellerStartEreduced(z,t1,t2,t3,t4,N): The expansion of the Even part of the #Feller generating function but z^2->z t1^2->t1 t3^2->t3 #For example, do: #FellerStartEreduced(z,t1,t2,t3,t4,4); FellerStartEreduced:=proc(z,t1,t2,t3,t4,N) local gu,i: gu:=Feller(z,t1,t2,t3,t4): gu:=normal(simplify((subs(z=-z,gu)+gu)/2)): gu:=expand(taylor(gu,z=0,2*N+6)): gu:=add(coeff(gu,z,2*i)*z^(i),i=0..N): gu:=add(coeff(gu,t1,2*i)*t1^(i),i=0..trunc((degree(gu,t1)+1)/2)): gu:=add(coeff(gu,t3,2*i)*t3^(i),i=0..trunc((degree(gu,t3)+1)/2)): gu:=expand(gu): end: