#Motzkin.txt
#This is the Motzkin analogue to Dr. Zeilberger's `Dyck.txt`
print(`To see the procedures in this package, type "Help()"`):
print(`For explanation of each procedure, type "WhatIs(procedure_name)"`):
print(`This package also contains procedures from Dr. Doron Zeilberger's FindRec.txt`):
print(`For a list of the procedures from FindRec.txt`):
print(` type "ezraFindRec();". For specific help type "ezraFindRec(procedure_name);" `):
print(`This package also contains procedures from Dr. Doron Zeilberger's SCHUTZENBERGER.txt`):
Help:=proc(): print(`SetU(m,n)`,`SetD(m,n)`,`SetF(m,n)`): print(`u(m,n)`,`d(m,n)`,`f(m,n)`):
print(`Ru(m,n,A,B,F)`, `Rd(m,n,A,B,F)`, `Rf(m,n,A,B,F)`):
print(`NRu(m,n,A,B,C,D,E,F)`, `NRd(m,n,A,B,C,D,E,F)`, `NRf(m,n,A,B,C,D,E,F)`):
print(`NRur(m,n,A,B,C,D,E,F,r)`, `NRdr(m,n,A,B,C,D,E,F,r)`, `NRfr(m,n,A,B,C,D,E,F,r)`):
print(`SeqABCDE(A,B,C,D,E,K)`, `SeqABCDEr(A,B,C,D,E,r,K)`):
print(`Motzkin(n)`, `Mk(n,k)`, `Targem1(L)`, `IsGoodCDE(L,C,D,E)`):
print(`Targem2(L)`, `IsGoodA(L,A)`, `IsGoodB(L,B)`, `MotzkinABCDE(A,B,C,D,E,n)`, `SeqABCDEbrute(A,B,C,D,E,K)`):
print(`IsGoodCDEr(L,C,D,E,r)`, `IsGoodAr(L,A,r)`, `IsGoodBr(L,A,r)`, `MotzkinPathsABCDEr(A,B,C,D,E,r,n)`, `SeqABCDErBrute(A,B,C,D,E,r,K)`):
print(`The following procedures from Dr. Doron Zeilberger's Dyck.txt are also included: IsRange1(f,k,n), IsRange(F,k,n), Aeq(gu,x,P)`):
end:
WhatIs:=proc(procedure):
if procedure=SetU then
print(`SetU(m,n): inputs non-negative integers m,n`):
print(`outputs set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step.`):
elif procedure=SetD then
print(`SetD(m,n): inputs non-negative integers m,n`):
print(`outputs set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step.`):
elif procedure=SetF then
print(`SetF(m,n): inputs non-negative integers m,n`):
print(`outputs set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step.`):
elif procedure=u then
print(`u(m,n): inputs non-negative integers m,n`):
print(`outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step.`):
elif procedure=d then
print(`d(m,n): inputs non-negative integers m,n`):
print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step.`):
elif procedure=f then
print(`f(m,n): inputs non-negative integers m,n`):
print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step.`):
elif procedure=Ru then
print(`Ru(m,n,A,B,F): inputs non-negative integers m and n, finite sets A and B of non-negative integers, and F={} or d or f`):
print(`Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step`):
print(`s.t. no peak heights belongs to A and no valley height belongs to B when F is appended to the walk.`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk.`):
elif procedure=Rd then
print(`Rd(m,n,A,B,F): inputs non-negative integers m and n, finite sets A and B of non-negative integers, and F={} or d or f`):
print(`Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step`):
print(`s.t. no peak heights belongs to A and no valley height belongs to B when F is appended to the walk.`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk.`):
elif procedure=Rf then
print(`Rf(m,n,A,B,F): inputs non-negative integers m and n, finite sets A and B of non-negative integers, and F={} or d or f`):
print(`Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step`):
print(`s.t. no peak heights belongs to A and no valley height belongs to B when F is appended to the walk.`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk.`):
elif procedure=NRu then
print(`NRu(m,n,A,B,C,D,E,F): inputs non-negative integers m and n, finite sets A, B,C,D, and E of non-negative integers, and F={} or d or f`):
print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step`):
print(`s.t. no peak has height belongs to A and no valley height belongs to B`):
print(`no upward run belongs to C, no downward run belongs to D, and no flat run belongs to E`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk.`):
elif procedure=NRd then
print(`NRd(m,n,A,B,C,D,E,F): inputs non-negative integers m and n, finite sets A, B,C,D, and E of non-negative integers, and F={} or u or f`):
print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step`):
print(`s.t. no peak has height belongs to A and no valley height belongs to B`):
print(`no upward run belongs to C, no downward run belongs to D, and no flat run belongs to E`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk.`):
elif procedure=NRf then
print(`NRf(m,n,A,B,C,D,E,F): inputs non-negative integers m and n, finite sets A, B,C,D, and E of non-negative integers, and F={} or u or d`):
print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step`):
print(`s.t. no peak has height belongs to A and no valley height belongs to B`):
print(`no upward run belongs to C, no downward run belongs to D, and no flat run belongs to E`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk.`):
elif procedure=IsRange1 then
print(`IsRange1(f,k,n): inputs a linear expression f in k and a non-negative integer n`):
print(`decides whether f(k0)=n for some non-negative integer k0.`):
print(`This was taken from Dr. Doron Zeilberger's Dyck.txt`):
elif procedure=IsRange then
print(`IsRange(F,k,n): inputs a set of linear expressions F={f} in k and a non-negative integer n`):
print(`decides whether f(k0)=n for some non-negative integer k0 and some f in F`):
print(`This was taken from Dr. Doron Zeilberger's Dyck.txt`):
elif procedure=NRur then
print(`NRur(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or d or f`):
print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step`):
print(`s.t. no peak has height is in the range of A and no valley height is in the range of B`):
print(`no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk.`):
elif procedure=NRdr then
print(`NRdr(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or u or f`):
print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step`):
print(`s.t. no peak has height is in the range of A and no valley height is in the range of B`):
print(`no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk.`):
elif procedure=NRfr then
print(`NRfr(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or u or d`):
print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step`):
print(`s.t. no peak has height is in the range of A and no valley height is in the range of B`):
print(`no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk.`):
elif procecedure=SeqABCDE then
print(`SeqABCDE(A,B,C,D,E,K): inputs finite sets A,B,C,D,E of non-negative integers and positive integer K`):
print(`Outputs the first K terms of the sequence of Motzkin paths of length n`):
print(`where no peak height belongs to A and no valley height belongs to B`):
print(`No upward run belongs to C, No downward run in D, and no flat run belongs to E`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk`):
elif procedure=SeqABCDEr then
print(`SeqABCDEr(A,B,C,D,E,r,K): inputs sets A,B,C,D,E of linear expressions in the variable r and positive integer K`):
print(`Outputs the first K terms of the sequence of Motzkin paths of length n`):
print(`where no peak has height is in the range of A and no valley height is in the range of B`):
print(`no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E`):
print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
print(`or a chain D F^k that occurs at the end of a walk`):
elif procedure=Motzkin then
print(`Motzkin(n): outputs all Motzkin paths of length n.`):
elif procedure=Targem1 then
print(`Targem1(L): translates a Motzkin path into a sequence of positive, negative integers, and sets`):
print(`representing the upwards, downwards, and flat runs respectively,`):
print(`such that the partial sums of the integers are always non-negative`):
elif procedure=IsGoodCDE then
print(`IsGoodCDE(L,C,D,E): Inputs a Motzkin path L, and finite sets C,D,E of non-negative integers.`):
print(`Decides whether no upwards run in L belongs to C, no downwards run in L belongs to D, and no flat run in L belongs to E.`):
elif procedure=IsGoodCDEr then
print(`IsGoodCDEr(L,C,D,E,r): Inputs a Motzkin path L, and sets C,D,E of linear expressions in the variable r.`):
print(`Decides whether no upwards run in L belongs to the range of C, no downward run in L belongs to the range of D, and no flat run in L belongs to the range of E.`):
elif procedure=Targem2 then
print(`Targem2(L): translates a Motzkin path into a sequence of the partial sums when it between a`):
print(`chain with steps in {0,1} starting with 1 (i.e. weakly increasing chain starting with an upwards step)`):
print(`and a chain with steps in {0,-1} starting with -1 (i.e. a weakly decreasing chain starting with a downwards step)`):
elif procedure=IsGoodA then
print(`IsGoodA(L,A): Inputs a Motzkin path L, and a finite set A of non-negative integers.`):
print(`Decides whether none of the heights in a peak location of L belongs to A`):
print(`A peak height is the maximal height of a chain beginning with the step 1 and ending with the step -1, s.t. any steps between (if any) are 0.`):
elif procedure=IsGoodB then
print(`IsGoodB(L,B): Inputs a Motzkin path L, and a finite set B of non-negative integers.`):
print(`Decides whether none of the heights in a valley location of L belongs to B.`):
print(`A valley height is the minimal height of a chain beginning with the step -1 and ending with the step 1 s.t. any steps between (if any) are 0.`):
elif procedure=IsGoodAr then
print(`IsGoodAr(L,A,r): Inputs a Motzkin path L, and a set A of linear expressions in the variable r.`):
print(`Decides whether none of the heights in a peak location belongs to the range of A`):
print(`A peak height is the maximal height of a chain beginning with the step 1 and ending with the step -1, s.t. any steps between (if any) are 0.`):
elif procedure=IsGoodBr then
print(`IsGoodBr(L,A,r): Inputs a Motzkin path L, and a set B of linear expressions in the variable r.`):
print(`Decides whether none of the heights in a valley location belongs to the range of B`):
print(`A valley height is the minimal height of a chain beginning with the step -1 and ending with the step 1 s.t. any steps between (if any) are 0.`):
elif procedure=MotzkinABCDE then
print(`MotzkinABCDE(A,B,C,D,E,n): Inputs five finite sets A, B, C,D,E of non-negative integers and a positive integer n.`):
print(`Outputs the set of Motzkin paths counted by SeqABCDE(A,B,C,D,E,K). In other words:`):
print(`where no peak can be of height that belongs to A and no valley can be of height that belongs to B`):
print(`No length of upward run belongs to C and no length of downward run belongs to D and no length of flat run belongs to E`):
print(`A peak height is the maximal height of a chain beginning with the step 1 (i.e. up) and ending with the step -1 (i.e. down), s.t. any steps between (if any) are 0 (i.e. flat).`):
elif procedure=MotzkinABCDEr then
print(`MotzkinABCDEr(A,B,C,D,E,r,n): Inputs five sets A, B, C,D,E of linear expressions in the variable r, and a positive integer n.`):
print(`Outputs the set of Motzkin paths counted by SeqABCDEr(A,B,C,D,E,r,K). In other words:`):
print(`where no peak can be of height that belongs to the range of A and no valley can be of height that belongs to the range of B,`):
print(`No length of upward run belongs to the range of C and no length of downward run belongs to the range of D and no length of flat run belongs to the range of E.`):
print(`A peak height is the maximal height of a chain beginning with the step 1 (i.e. up) and ending with the step -1 (i.e. down), s.t. any steps between (if any) are 0 (i.e. flat).`):
elif procedure=SeqABCDEbrute then
print(`SeqABCDEbrute(A,B,C,D,E,K): Same output as SeqABCDE(A,B,C,D,E,K) but by complete brute force.`):
print(`FOR CHECKING PURPOSES ONLY; Don't make K too big`):
elif procedure=SeqABCDErBrute then
print(`SeqABCDErBrute(A,B,C,D,K): Same output as SeqABCDEr(A,B,C,D,E,r,K) but by complete brute force.`):
print(`FOR CHECKING PURPOSES ONLY; Don't make K too big.`):
elif procedure=Aeq then
print(`Aeq(L,x,P): The guessed algebraic equation in P(x) satisfied by the generating function of the sequence whose beginning`):
print(`starting at n=0 is L. Try:`):
print(`Aeq([seq(binomial(2*i,i)/(i+1),i=1..10)],x,P);`):
print(`This was taken from Dr. Doron Zeilberger's Dyck.txt`):
fi:
end:
#SetU(m,n): inputs non-negative integers m,n
#Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step.
SetU:=proc(m,n) local gu,mu,r,mu1 :
option remember:
if (m<0 or n<0 or m=0 and type(k0,integer) then
true:
else
false:
fi:
end:
#IsRange(F,k,n): inputs a set of linear expressions F={f} in k and a non-negative integer n
#decides whether f(k0)=n for some non-negative integer k0 and some f in F
#This was taken from Dr. Doron Zeilberger's Dyck.txt
IsRange:=proc(F,k,n) local f:
for f in F do
if IsRange1(f,k,n) then
RETURN(true):
fi:
od:
false:
end:
#NRur(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or d or f
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step
#s.t. no peak has height is in the range of A and no valley height is in the range of B
#no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk
NRur:=proc(m,n,A,B,C,D,E,F,r) local gu,k:
option remember:
if (m<0 or n<0 or m1 then
if type(L1[i-1],positive)=type(L1[i+1],positive) then
su:=su+L1[i+1]:
L2:=[op(1..-2,L2),su]:
i:=i+1:
fi: fi: fi: od:
L2:
end:
#IsGoodA(L,A): Inputs a Motzkin path L, and a finite set A of non-negative integers.
#Decides whether none of the heights in a peak location of L belongs to A
#A peak height is the maximal height of a chain beginning with the step 1 and ending with the step -1,
#s.t. any steps between (if any) are 0.
IsGoodA:=proc(L,A) local L1,i:
L1:=Targem2(L):
for i from 1 to nops(L1)/2 do
if member(L1[2*i-1],A) then
RETURN(false):
fi:
od:
true:
end:
#IsGoodB(L,B): Inputs a Motzkin path L, and a finite set B of non-negative integers,
#Decides whether none of the heights in a valley location of L belongs to B.
#A valley height is the minimal height of a chain beginning with the step -1 and ending with the step 1,
#s.t. any steps between (if any) are 0.
IsGoodB:=proc(L,B) local L1,i:
L1:=Targem2(L):
for i from 1 to nops(L1)/2 do
if member(L1[2*i],B) then RETURN(false): fi:
od:
true:
end:
#IsGoodAr(L,A,r): Inputs a Motzkin path L, and a set A of linear expressions in the variable r.
#Decides whether none of the heights in a peak location belongs to the range of A
#A peak height is the maximal height of a chain beginning with the step 1 and ending with the step -1,
#s.t. any steps between (if any) are 0.
IsGoodAr:=proc(L,A,r) local L1,i:
L1:=Targem2(L):
for i from 1 to nops(L1)/2 do
if IsRange(A,r,L1[2*i-1]) then
RETURN(false):
fi:
od:
true:
end:
#IsGoodBr(L,A,r): Inputs a Motzkin path L, and a set B of linear expressions in the variable r.
#Decides whether none of the heights in a valley location belongs to the range of B
#A valley height is the minimal height of a chain beginning with the step -1 and ending with the step 1,
#s.t. any steps between (if any) are 0.
IsGoodBr:=proc(L,B,r) local L1,i:
L1:=Targem2(L):
for i from 1 to nops(L1)/2 do
if IsRange(B,r,L1[2*i]) then
RETURN(false):
fi:
od:
true:
end:
#MotzkinABCDE(A,B,C,D,E,n): Inputs five finite sets A, B, C,D,E of non-negative integers and a positive integer n,
#Outputs the set of Motzkin paths counted by SeqABCDE(A,B,C,D,E,K). In other words:
#where no peak can be of height that belongs to A and no valley can be of height that belongs to B
#No length of upward run belongs to C and no length of downward run belongs to D and no length of flat run belongs to E
#A peak height is the maximal height of a chain beginning with the step 1 (i.e. up) and ending with the step -1 (i.e. down),
#s.t. any steps between (if any) are 0 (i.e. flat).
MotzkinABCDE:=proc(A,B,C,D,E,n) local gu,mu,mu1:
mu:=Motzkin(n):
gu:={}:
for mu1 in mu do
if (IsGoodCDE(mu1,C,D,E) and IsGoodA(mu1,A) and IsGoodB(mu1,B)) then
gu:=gu union {mu1}:
fi:
od:
gu:
end:
#MotzkinABCDEr(A,B,C,D,E,r,n): Inputs five sets A, B, C,D,E of linear expressions in the variable r, and a positive integer n.
#Outputs the set of Motzkin paths counted by SeqABCDEr(A,B,C,D,E,r, K). In other words:
#where no peak can be of height that belongs to the range of A and no valley can be of height that belongs to the range of B,
#No length of upward run belongs to the range of C and no length of downward run belongs to the range of D
#and no length of flat run belongs to the range of E.
#A peak height is the maximal height of a chain beginning with the step 1 (i.e. up) and ending with the step -1 (i.e. down),
#s.t. any steps between (if any) are 0 (i.e. flat).
MotzkinABCDEr:=proc(A,B,C,D,E,r,n) local gu,mu,mu1:
mu:=Motzkin(n):
gu:={}:
for mu1 in mu do
if (IsGoodCDEr(mu1,C,D,E,r) and IsGoodAr(mu1,A,r) and IsGoodBr(mu1,B,r)) then
gu:=gu union {mu1}:
fi:
od:
gu:
end:
#SeqABCDEbrute(A,B,C,D,E,K): Same output as SeqABCDE(A,B,C,D,E,K) but by complete brute force.
#FOR CHECKING PURPOSES ONLY
#Don't make K too big.
SeqABCDEbrute:=proc(A,B,C,D,E,K) local i:
[seq(nops(MotzkinABCDE(A,B,C,D,E,i)),i=0..K)]:
end:
#SeqABCDErBrute(A,B,C,D,K): Same output as SeqABCDEr(A,B,C,D,E,r,K) but by complete brute force.
#FOR CHECKING PURPOSES ONLY
#Don't make K too big.
SeqABCDErBrute:=proc(A,B,C,D,E,r,K) local i:
[seq(nops(MotzkinABCDEr(A,B,C,D,E,r,i)),i=0..K)]:
end:
###Start from SCHUTZENBERGER.txt
#empir(gu,degx,degP,x,P)
#to "fit" an algebraic equation F(P(x),x)=0 of degree
#degP in P(x) and degx in n for P(x):=sum_i gu[i]*x^i
empir:=proc(gu,degx,degP,x,P)
local i1,i2,F,a,cand,lu,eq,var,mu,flo,pip,ka,mu1,Ftry:
if (1+degx)*(1+degP) > nops(gu)-15 then
RETURN(`sequence too small`):
fi:
F:=0:
var:={}:
for i1 from 0 to degx do
for i2 from 0 to degP do
F:=F+a[i1,i2]*x^i1*P^i2:
var:=var union {a[i1,i2]}:
od:
od:
cand:=0:
for i1 from 0 to nops(gu)-1 do
cand:=cand+op(i1+1,gu)*x^i1:
od:
lu:=subs(P=cand,F):
lu:=taylor(lu,x=0,nops(gu)-1):
eq:={}:
for i1 from 0 to nops(gu)-2 do
eq:=eq union {coeff(lu,x,i1)=0}
od:
mu:=solve(eq,var):
ka:=0:
for mu1 in mu do
if op(1,mu1)=op(2,mu1) then
ka:=ka+1:
fi:
od:
if ka>1 then
RETURN(FAIL):
fi:
F:=subs(mu,F):
if F=0 then
RETURN(FAIL):
fi:
flo:=degree(F,P):
pip:=coeff(F,P,flo):
flo:=degree(pip,x):
pip:=coeff(pip,x,flo):
F:=normal(F/pip):
Ftry:=taylor(subs(P=add(gu[i1+1]*x^i1,i1=0..nops(gu)-1),F),x=0,nops(gu)+1):
if {seq(coeff(Ftry,x,i1),i1=0..nops(gu)-2)}<>{0} then
RETURN(FAIL):
fi:
normal(F):
end:
Aeq:=proc(gu,x,P)
local degx,degP,L,lu:
for L from 1 to (nops(gu)-15)/3 do
for degP from 1 to L do
for degx from 0 to min(trunc((nops(gu)-15)/(1+degP))-1,L-degP) do
lu:=empir(gu,degx,degP,x,P):
if lu<>FAIL then
RETURN(lu):
fi:
od:
od:
od:
FAIL:
end:
###end from SCHUTZENBERGER.txt
##start from Findrec
ezraFindrec:=proc()
if args=NULL then
print(` FindRec.txt: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures: `):
print(` findrec, Findrec, FindrecF`):
print():
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):
RETURN(FAIL):
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):
RETURN(FAIL):
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
ope:=ope[1]:
if {seq(add(subs(n=n1+i1,coeff(ope,N,i1))*f[n1+i1],i1=0..degree(ope,N)),n1=1..nops(f)-degree(ope,N))}<>{0} then
RETURN(FAIL):
fi:
ope:=Yafe(ope,N)[2]:
if SeqFromRec(ope,n,N,[op(1..degree(ope,N),f)],nops(f))<>f then
RETURN(FAIL):
else
RETURN(ope):
fi:
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 and degree(ope,N)<>0 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,ka,i1:
ope1:=ope:
L:=degree(ope1,N):
if nops(Ini)<>L then
ERROR(`Ini should be of length`, L):
fi:
ka:={seq(subs(n=i1,lcoeff(ope1,N)),i1=1..K)}:
if member(0,ka) then
RETURN(FAIL):
fi:
ope1:=expand(subs(n=n-L,ope1)/N^L):
gu:=Ini:
for n1 from nops(Ini)+1 to K do
if member(0, { seq( subs(n=n1,denom(coeff(ope1,N,-j1) )) ,j1=1..L )}) then
RETURN(FAIL):
else
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
fi:
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):
RETURN(FAIL):
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 E10 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