###################################################################### ## CTcong.txt Save this file as CTcong.txt stay in the same # ## directory, get into Maple (by typing: maple ) # ## and then type: read CTcong.txt # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## zeilberg at math dot rutgers dot edu. # ###################################################################### print(`This is CTcong.txt, A Maple package`): print(`accompanying an article by William Y.C. Chen, Qing-Hu Hou and Doron Zeilberger entitled: `): print(`"Automated Discovery and Proof of Congruence Theorems for Partial Sums of Combinatorial Sequences" `): print(): print(`The most current version is available on WWW at:`): print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/CTcong.txt`): print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`): print(): print(`-----------------------------------`): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): print(`-----------------------------------`): print(`For a list of the supporting functions type: ezra1();`): print(` For specific help type "ezra(procedure_name);" `): print(): print(`-----------------------------------`): print(`-----------------------------------`): ezraOne:=proc() if args=NULL then print(`The old, deperated, procedures, are: CTtoSeq, CTtoSeqModp, EvSum, EvSumC, SeqCT, Theo `): print(``): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(`The SUPPORTING procedures are: EvalQP, CHZcheck , EvSumCheck, PaSu1check, EvSumSep, EvSumSepG, RatToQP, SeqC, S1 `): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` CTcong: A Maple package for automatically discovering and proving congruence theorem for partial sums of combinatorial sequences given by Constant Terms Expressions.`): print(`The MAIN procedures are: SEQ, PaSu, PaSu1, CHZ, TheoG, TheoQP `): print(``): elif nargs=1 and args[1]=CTtoSeq then print(`CTtoSeq(P,Q,x,N): The first N terms (starting at 0) of the sequence CT(Q(x)*P(x)^n)), where Q and P Are Laurent polynomials in x. Try`): print(`CTtoSeq((1+x)^2/x,1,x,10);`): elif nargs=1 and args[1]=CTtoSeqModp then print(`CTtoSeqModp(P,Q,x,N,p): The first N terms (starting at 0) of the sequence CT(Q(x)*P(x)^n)) mod p, where Q and P Are Laurent polynomials in x, and p is a prime. Try`): print(`CTtoSeqModp((1+x)^2/x,1,x,10,101);`): elif nargs=1 and args[1]=EvalQP then print(`EvalQP(Q,n,n0): evaluates the quasi-polynomial Q of n at n0`): print(` For example try: `): print(` EvalQP([2*n^3,n^3],n,5); `): elif nargs=1 and args[1]=EvSum then print(`EvSum(P,Q,x,p): inputs Laurent polynomials P,Q, (with integer coefficients)in the variable x (and the coefficient of the lowest degree term of P must be 1 )`): print(` and a symbol p (denoting a general prime), outputs`): print(`a simplified expression for SeqCT(P,Q,x,N) `): print(`For example, try:`): print(`EvSum(1/x+2+x,1,x,p); `): print(`EvSum(1/x+1+x,1-x^2,x,p); `): print(`EvSum(1/x+3+2*x,1,x,p); `): elif nargs=1 and args[1]=EvSumC then print(`EvSumC(P,Q,x,p,C): Like EvSum(P,Q,x,p) but also inputs a symbol, C, where C(i) denotes the coefficient of the outputted rational function (the second part of the output)`): print(`the first part in a linear compination of terms of the form C(i*p-j, for specific integers i and j. Try:`): print(`EvSumC((1/x+2+x),1,x,p,C);`): print(`EvSumC(1/x+2+x,1-x,x,p,C); `): print(`EvSumC(1/x+1+x,1-x^2,x,p,C); `): print(`EvSumC(1/x+3+2*x,1,x,p,C); `): elif nargs=1 and args[1]=PaSu then print(`PaSu(P,Q,x,p,C,r): Inputs Laurent polynomials P,Q, (with integer coefficients)in the variable x (and the coefficient of the lowest degree term of P must be 1 ) `): print(` and a symbol p (denoting a general prime), as well as a letter C, and a positive integer r. Outputs`): print(`a pair whose second component is a rational function of x, and denoting by C(i) the coefficient of x^i in its Maclaurin expansion, the first component is`): print(` an expression that equals, mod p, to the partial sum from 0 to r*p-1, of the combinatorial sequence given by`): print(`A(i):=ConstandTermOf Q(x)*P(x)^i `): print(`PaSu((1/x+2+x),1,x,p,C,2);`): print(`PaSu(1/x+2+x,1-x,x,p,C,3); `): print(`PaSu(1/x+1+x,1-x^2,x,p,C,4); `): print(`PaSu(1/x+3+2*x,1,x,p,C,2); `): elif nargs=1 and args[1]=PaSu1 then print(`PaSu1(P,Q,x,p,r): inputs Laurent polynomials P,Q, (with integer coefficients)in the variable x (and the coefficient of the lowest degree term of P must be 1) `): print(` and a symbol p (denoting a general prime), as well as a positive integer r, outputs`): print(` an expression whose constant term, mod p, is the partial sum from 0 to r*p-1, mod p, of the combinatorial sequence given by`): print(`A(i):=ConstandTermOf Q(x)*P(x)^i `): print(`For example, try:`): print(`PaSu1(1/x+2+x,1,x,p,2); `): print(`PaSu1(1/x+1+x,1-x^2,x,p,3); `): print(`PaSu1(1/x+3+2*x,1,x,p,4); `): elif nargs=1 and args[1]=PaSu1check then print(`PaSu1check(P,Q,x,N,r): checks the output for PaSu1(P,Q,x,p,r) (or, equivalently, PaSu(P,Q,x,p,C,r), up to all primes<=N.`): print(` Try: `): print(` PaSu1check(1/x+2+x,1,x,200,2); `): elif nargs=1 and args[1]=EvSumCheck then print(`EvSumCheck(P,Q,x,N): checks the output for EvSum(P,Q,x,p) (or, equivalently, EvSumC(P,Q,x,p,C), up to all primes<=N.`): print(` Try: `): print(` EvSumCheck(1/x+2+x,1,x,200); `): elif nargs=1 and args[1]=EvSumSep then print(`EvSumSep(P,Q,x,y): Like EvSum(P,Q,x,p) (q.v.), but with x^p, renamed 1/y, part separated. Returns the pair [(x^p stuff), Ratioanl funcion]`): print(`For example, try:`): print(`EvSumSep(1/x+2+x,1,x,y); `): print(`EvSumSep(1/x+1+x,1-x^2,x,y); `): print(`EvSumSep(1/x+3+2*x,1,x,y); `): elif nargs=1 and args[1]=EvSumSepG then print(`EvSumSepG(P,Q,x,y,r): Like PaSu1(P,Q,x,p,r) (q.v.), but with x^p, renamed 1/y, part separated. Returns the pair [(x^p stuff), Ratioanl funcion]`): print(`For example, try:`): print(`EvSumSepG(1/x+2+x,1,x,y,2); `): print(`EvSumSepG(1/x+1+x,1-x^2,x,y,3); `): print(`EvSumSep(1/x+3+2*x,1,x,y,4); `): elif nargs=1 and args[1]=EvQPg then print(`EvQPg(P,Q,x,p,C,r,d,n): Like PaSu(P,Q,x,p,C,r), but expresses C)i) as A quasi-polynomial of degree d, in n, rather than a sequence of coefficients of a rational function.`): print(`Try:`): print(` EvQPg((1/x+2+x),1,x,p,C,2,0,n); `): print(`EvQPg(1/x+1+x,1-x^2,x,p,C,3,0,n); `): print(`EvQPg(1/x^2+1/x+1+x+x^2,1-x^2,x,p,C,3,0,n); `): elif nargs=1 and args[1]=CHZ then print(`CHZ(P,Q,x,r,d,p): The Chen-Hou-Zeilberger algorithm. Inputs Laurent polynomials, P and Q, with integer coefficients, of the variable x, as well as a positive integer r`): print(`another positive integer d, and a symbol p, denoting a general (symbolic) prime, and outputs, if found, a quasi-polynomial expression for the expression`): print(`Sum(A(i),i=0..r*p-1) mod p `): print(` where A(i) is the combinatorial sequence defined by the Constant Term Expression `): print(`A(i):=ConstandTermOf Q(x)*P(x)^i `): print(`The quasi-polynomial is expressed as a list of expressions in p (of course sometimes they may be constant), meaning that value is given by the i-th component if the prime is`): print(`i modulo the length of the list. `): print(`For example, [p^2,p,4*p-2] represents the quasi-polynmial of quasi-period 3 that equals p^2 if p is 1 mod 3, equals p, if p is 2 mod 3, and equals 4*p-2 if p is 0 mod 3`): print(` For the sum of the Central Binomial Coefficients from 0 to r*p-1 for r=1,..., 5, type, in order`): print(` CHZ(1/x+2+x,1,x,1,0,p); `): print(` CHZ(1/x+2+x,1,x,2,0,p); `): print(` CHZ(1/x+2+x,1,x,3,0,p); `): print(` CHZ(1/x+2+x,1,x,4,0,p); `): print(` CHZ(1/x+2+x,1,x,5,0,p); `): print(` For the sum of the Catalan Numbers from 0 to r*p-1 for r=1,..., 5, type, in order`): print(` CHZ(1/x+2+x,1-x,x,1,0,p); `): print(` CHZ(1/x+2+x,1-x,x,2,0,p); `): print(` CHZ(1/x+2+x,1-x,x,3,0,p); `): print(` CHZ(1/x+2+x,1-x,x,4,0,p); `): print(` CHZ(1/x+2+x,1-x,x,5,0,p); `): print(` For the sum of the Motzkin Numbers from 0 to r*p-1 for r=1,..., 5, type, in order`): print(` CHZ(1/x+1+x,1-x^2,x,1,0,p); `): print(` CHZ(1/x+1+x,1-x^2,x,2,0,p); `): print(` CHZ(1/x+1+x,1-x^2,x,3,0,p); `): print(` CHZ(1/x+1+x,1-x^2,x,4,0,p); `): print(` CHZ(1/x+1+x,1-x^2,x,5,0,p); `): print(` For the sum of the Central Pentagonal Numbers, from 0 to r*p-1 for r=1,..., 5, type, in order`): print(` CHZ(1/x^2+1/x+1+x+x^2,1,x,1,1,p); `): print(` CHZ(1/x^2+1/x+1+x+x^2,1,x,2,1,p); `): print(` CHZ(1/x^2+1/x+1+x+x^2,1,x,3,1,p); `): print(` CHZ(1/x^2+1/x+1+x+x^2,1,x,4,1,p); `): print(` CHZ(1/x^2+1/x+1+x+x^2,1,x,5,1,p); `): elif nargs=1 and args[1]=CHZcheck then print(`CHZcheck(P,Q,x,r,d,N): Nunerically checks CHZ(P,Q,x,r,d,p) for all primes <=N.`): print(`Try: `): print(` CHZcheck((1/x+2+x),1,x,3,0,200); `): print(` CHZcheck((1/x+1+x),1-x^2,x,5,0,100); `): elif nargs=1 and args[1]=RatToQP then print(`RatToQP(R,x,L,d,n): inputs a rational function R, in the variable x, a positive integer L, and a degree d, tries to find`): print(`a quasi-polynomial of degree d in n, describing its coefficients. Try:`): print(`RatToQP(1/(1+x^3)^2,x,100,1,n);`): elif nargs=1 and args[1]=S1 then print(` S1(F,k,p): inputs a discrete expression in k, F, that evaluates to integers and finds the partial sum to p-1 mod p. Try: `): print(` S1(binomial(2*k,k),k,11); `): elif nargs=1 and args[1]=SeqC then print(`SeqC(F,k,N): Inputs a hypergeoemetric term the first N terms of the sequence`): print(` Sum(F(k),k=0..p-1) mod p for the first N primes (starting at 5). Try: `): print(` SeqC(binomial(2*k,k),k,20): `): elif nargs=1 and args[1]=SeqCT then print(`SeqCT(P,Q,x,N): Inputs Laurent polynomials P and Q in the variable x and a positive integer N, outputs `): print(`Sum(F(k),k=0..p-1) mod p for the first N primes (starting at 5) `): print(`where F(k) is the Constant Term of Q(x)*P(x)^k`): print(`SeqCT((1+x)^2/x,1,x,40); `): print(`SeqCT((1+x)^2/x,1-x,x,40); `): print(`SeqCT(1/x+1+x,1-x^2,x,40); `): print(`SeqCT(1/x+3+2*x,1,x,40); `): elif nargs=1 and args[1]=SEQ then print(`SEQ(P,Q,x,N,r): Inputs Laurent polynomials P and Q in the variable x and a positive integer N, as well as a positive integer r; outputs `): print(`Sum(F(k),k=0..r*p-1) mod p (where moduli larger than p/2 are denoted by negative integers (e.g. 4 mod 5 is -1)), for the first N primes (starting at 5) `): print(`where F(k) is the Constant Term of Q(x)*P(x)^k`): print(`SEQ((1+x)^2/x,1,x,40,3); `): print(`SEQ((1+x)^2/x,1-x,x,40,2); `): print(`SEQ(1/x+1+x,1-x^2,x,40,4); `): print(`SEQ(1/x+3+2*x,1,x,20,5); `): elif nargs=1 and args[1]=Theo then print(`Theo(P,Q,x,p,C): inputs Laurent polynomials P and Q, in the variable x, and symbols C and p, outputs a theorem`): print(`about the partial sum (up to p-1) mod p of the combinatorial sequence defined by C(n):=ConstantTerm (Q*P^n).`): print(`Try: `): print(`Theo(1/x+2+x,1,x,p,C):`): elif nargs=1 and args[1]=TheoG then print(`TheoG(P,Q,x,p,C,r): inputs Laurent polynomials P and Q, in the variable x, and symbols C and p, as well as a positive integer r, outputs a theorem`): print(`about the partial sum (up to r*p-1) mod p of the combinatorial sequence defined by C(n):=ConstantTerm (Q*P^n).`): print(`Try: `): print(`TheoG(1/x+2+x,1,x,p,C,2):`): elif nargs=1 and args[1]=TheoQP then print(`TheoQP(P,Q,x,p,r,d): inputs Laurent polynomials P and Q, in the variable x, and symbols C and p, and a positive integer r,outputs a theorem`): print(`about the partial sum (up to p-1) mod p of the combinatorial sequence defined by C(n):=ConstantTerm (Q*P^n).`): print(`expressing it in terms of quasipolynomials of degree<=d `): print(`Of course, it may not be possible, then it returns FAIL. TO get a guaranteed theorem use procedure TheoG `): print(`Try: `): print(`For the central binomial coefficients from 0 to 2p-1 try: `): print(`TheoQP(1/x+2+x,1,x,p,2,0):`): print(`Try: `): print(`For the Catalan numbers from 0 to 3p-1 try: `): print(`TheoQP(1/x+2+x,1-x,x,p,3,0):`): print(`For the Motzkin numbers from 0 to 5p-1 try: `): print(`TheoQP(1/x+1+x,1-x^2,x,p,5,0):`): print(`For the Central Pentagonal numbers from 0 to p-1 try: `): print(`TheoQP(1/x^2+1/x+1+x+x^2,1,x,p,1,1):`): print(`For the Central Septagonal numbers from 0 to p-1 try: `): print(`TheoQP(1/x^3+1/x^2+1/x+1+x+x^2+x^3,1,x,p,2,0):`): else print(`There is no such thing as`, args): fi: end: ######stat Guesing Quasi-polynomials ##begin guessing polynomials #GuessPol1(L,d,n): guesses a polynomial of degree d in n for # the list L, such that P(i)=L[i+1] for i=1..nops(L) #For example, try: #GuessPol1([seq(i,i=1..10)],1,n); GuessPol1:=proc(L,d,n) local P,i,a,eq,var: if d>nops(L)-2 then ERROR(`the list is too small`): fi: P:=add(a[i]*n^i,i=0..d): var:={seq(a[i],i=0..d)}: eq:={seq(subs(n=i,P)-L[i],i=1..nops(L))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: subs(var,P): end: #GuessPol(L,n): guesses a polynomial of degree d in n for # the list L, such that P(i)=L[i] for i=1..nops(L) for example, try: #GuessPol([seq(i,i=1..10)],n); GuessPol:=proc(L,n) local d,gu: for d from 0 to nops(L)-2 do gu:=GuessPol1(L,d,n): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: ##end guessing polynomials #GuessQP1(L,r,n): Inputs a list L and outputs a list of polynomials #of length r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for #i=1..r, For example try: #GuessQP1:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],2,n); GuessQP1:=proc(L,r,n) local M,i,L1,gu,m: M:=[]: for i from 1 to r do L1:=[seq(L[i+(m-1)*r],m=1..trunc((nops(L)-i)/r)+1)]: gu:=GuessPol(L1,m): if gu=FAIL then RETURN(FAIL): fi: gu:=expand(subs(m=(n-i)/r+1,gu)): M:=[op(M),gu]: end: M: end: #GuessQP11(L,r,n,d): Inputs a list L, positive integer r, symbol n, and a positive integer dand outputs a list of polynomials #of degree r length r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for #i=1..r, For example try: #GuessQP11:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],2,n,0); GuessQP11:=proc(L,r,n,d) local M,i,L1,gu,m: M:=[]: for i from 1 to r do L1:=[seq(L[i+(m-1)*r],m=1..trunc((nops(L)-i)/r)+1)]: gu:=GuessPol1(L1,d,m): if gu=FAIL then RETURN(FAIL): fi: gu:=expand(subs(m=(n-i)/r+1,gu)): M:=[op(M),gu]: end: M: end: #GuessQP(L,n): Inputs a list L and a variable n, outputs a list of polynomials #of length r, for some r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for #i=1..r, For example try: #GuessQP:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],n); GuessQP:=proc(L,n) local r,gu: for r from 1 to trunc(nops(L)/4) do gu:=GuessQP1(L,r,n): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: #GuessQPd(L,n,d): Inputs a list L and a variable n, and positive integer d,outputs a list of degree-d polynomials #of length r, for some r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for #i=1..r, For example try: #GuessQPd:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],n,1); GuessQPd:=proc(L,n,d) local r,gu: for r from 1 to trunc(nops(L)/4) do gu:=GuessQP11(L,r,n,d): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: #EvalQP(Q,n,n0): evaluates the quasi-polynomial Q of n at n0 #For example try: #EvalQP([2*n^3,n^3],n,5); EvalQP:=proc(Q,n,n0) local i: i:=n0 mod nops(Q): if i=0 then i:=nops(Q): fi: subs(n=n0,Q[i]): end: #RatToQP(R,x,L,d,n): inputs a rational function R, in the variable x, a positive integer L, and a degree d, tries to find #a quasi-polynomial describing its coefficients. Try: #RatToQP(1/(1+x^3)^2,x,100,1); RatToQP:=proc(R,x,L,d,n) local lu,R1: R1:=rem(numer(R),denom(R),x)/denom(R): lu:=taylor(R1,x=0,L+1): lu:=[seq(coeff(lu,x,i),i=1..L)]: GuessQPd(lu,n,d): end: ######end Guesing Quasi-polynomials #S1(F,k,p): inputs a discrete expression in k, F, that evaluates to integers and finds the partial sum to p-1 mod p. Try: #S1(binomial(2*k,k),k,11); S1:=proc(F,k,p) local i: mods(add(eval(subs(k=i,F)),i=0..p-1),p): end: #SeqC(F,k,N): Inputs a hypergeoemetric term the first N terms of the sequence #Sum(F(k),k=0..p-1) mod p for the first N primes (starting at 5) . Try: #SeqC(binomial(2*k,k),k,p,N): SeqC:=proc(F,k,N) local i,gu: gu:=[seq( [ithprime(i), S1(F,k,ithprime(i))],i=3..N)]: gu,{seq(gu[i][2],i=1..nops(gu))} : end: #S1c(F,k,p): inputs a discrete expression in k, F, that evaluates to integers and finds the partial sum to p-1 mod p. Try: #S1(binomial(2*k,k),k,11); S1c:=proc(F,k,p) local i: mods(add(eval(subs(k=i,F)),i=0..p-1),p): end: #CTtoSeq:(P,Q,x,N): The first N terms (starting at 0) of the sequence CT(Q(x)*P(x)^n)), where Q and P Are Laurent polynomials in x. Try #CToSeq((1+x)^2,1,x,10); CTtoSeq:=proc(P,Q,x,N) local gu,L,i: gu:=Q: L:=[coeff(gu,x,0)]: for i from 1 to N do gu:=expand(gu*P): L:=[op(L),coeff(gu,x,0)]: od: L: end: #SeqCT(P,Q,x,N): Inputs Laurent polynomials P and Q in the variable x and a positive integer N, outputs #Sum(F(k),k=0..p-1) mod p for all primes less than N #where F(k) is the Constant Term of Q(x)*P(x)^k. Try: #SeqCT((1/x+1+x,1,x,20); SeqCT:=proc(P,Q,x,N) local L,p,gu,mu,i: L:=CTtoSeq(P,Q,x,N+1): p:=5: gu:=[]: while p<=N do mu:=mods(convert([op(1..p,L)],`+`),p): gu:=[op(gu),[p,mu]]; p:=nextprime(p): od: gu,{seq(gu[i][2],i=1..nops(gu))} : end: #EvSum(P,Q,x,p): inputs Laurent polynomials P,Q, (with integer coefficients)in the variable x, and a symbol p (denoting a general prime), outputs #a simplified expression for SeqCT(P,Q,x,N) EvSum:=proc(P,Q,x,p) local m,POL,POL1,i: m:=-ldegree(P,x): if coeff(P,x,-m)<>1 then print(`The coefficient of the lowest power of`, P, `should have been 1 `): RETURN(FAIL): fi: POL:=expand(x^m*P): POL1:=add(coeff(POL,x,i)*x^i, i=0..m-1): POL1:=subs(x=x^p,POL1): Q*POL1/(POL-x^m)/x^(m*(p-1)): end: #EvSumSep(P,Q,x,y): Like EvSum(P,Q,x,p) but with the rational function part seperated. Outputs it in the form #where y stands for 1/x^p: #[Expression with x^p, RationalFunction of x]. EvSumSep:=proc(P,Q,x,y) local m,POL,POL1,i: option remember: m:=-ldegree(P,x): if coeff(P,x,-m)<>1 then print(`The coefficient of the lowest power of`, P, `should have been 1 `): RETURN(FAIL): fi: POL:=expand(x^m*P): POL1:=add(coeff(POL,x,i)*x^i, i=0..m-1): POL1:=subs(x=1/y,POL1): #Q*POL1/(POL-x^m)/x^(m*(p-1)): [expand(POL1*x^m*y^m), factor(Q/(POL-x^m))] end: #EvSumC(P,Q,x,p,C): Like EvSum(P,Q,x,p) but also inputs a symbol, C, where C(i) denotes the coefficient of the outputted rational function (the second part of the output) #the first part in a linear compination of terms of the form C(i*p-j, for specific integers i and j. Try: #EvSumC((1/x+2+x),1,x,p,C); EvSumC:=proc(P,Q,x,p,C) local y,gu,mu,R,gu1,i,j: gu:=EvSumSep(P,Q,x,y): if gu=FAIL then RETURN(FAIL): fi: R:=gu[2]: gu:=gu[1]: mu:=0: for i from 0 to degree(gu,y) do gu1:=coeff(gu,y,i): for j from ldegree(gu1,x) to degree(gu1,x) do mu:=mu+C(i*p-j)*coeff(gu1,x,j): od: od: [mu,R]: end: #Theo(P,Q,x,p,C): inputs Laurent polynomials P and Q, in the variable x, and symbols C and p, outputs a theorem #about the partial sum (up to p-1) mod p of the combinatorial sequence defined by C(n):=ConstantTerm (Q*P^n). #Try: #Theo(1/x+2+x,1,x,p,C): Theo:=proc(P,Q,x,p,C) local gu,A,i: gu:=EvSumC(P,Q,x,p,C): if gu=FAIL then RETURN(FAIL): fi: print(``): print(`--------------------------------------------------------`): print(``): print(`Theorem: Define the following sequence: A(i) is the Constant term of the Laurent polynomial `): print(Q*P^i): print(`and for any prime p, let`): print(B(p)=Sum(A(i),i=0..p-1)): print(`Let the C-finite sequence C(n) be defined in terms of the generating function`): print(``): print(Sum(C(i)*x^i,i=0..infinity)=gu[2]): print(`then, modulo p, `): print(B(p)=gu[1]): print(``): end: #CTtoSeqModp(P,Q,x,N,p): The first N terms (starting at 0) of the sequence CT(Q(x)*P(x)^n)) mod p, where Q and P Are Laurent polynomials in x, and #p is a prime (in fact any integer). Try #CToSeq(Modp(1+x)^2,1,x,10,101); CTtoSeqModp:=proc(P,Q,x,N,p) local gu,L,i: gu:=Q: L:=[coeff(gu,x,0)] mod p: for i from 1 to N do gu:=expand(gu*P) mod p: L:=[op(L),coeff(gu,x,0)]: od: L: end: #EvSumCheck(P,Q,x,N): checks the output for EvSum(P,Q,x,p) (or, equivalently, EvSumC(P,Q,x,p,C), up to all primes<=N. #Try: #EvSumCheck(1/x+2+x,1,x,200); EvSumCheck:=proc(P,Q,x,N) local y,gu,lu,mu,i,j,mu1,p,coe1: gu:=EvSumSep(P,Q,x,y): lu:=taylor(gu[2],x=0,N*degree(gu[1],y)+1): lu:=[seq(coeff(lu,x,i),i=1..N*degree(gu[1],y))]: mu:=[]: p:=5: while p<=N do mu1:=0: for i from 1 to degree(gu[1],y) do for j from 0 to degree(coeff(gu[1],y,i),x) do coe1:=coeff(coeff(gu[1],y,i),x,j): mu1:=mu1+coe1*lu[i*p-j]: od: od: mu:=[op(mu),[p, mods(mu1,p)]]: p:=nextprime(p): od: evalb(mu=SeqCT(P,Q,x,N)[1]): end: ##srart generalized procedured #SEQ(P,Q,x,N,r): Inputs Laurent polynomials P and Q in the variable x and a positive integer N, and a positive integer r, outputs #Sum(F(k),k=0..r*p-1) mod p for all primes less than N #where F(k) is the Constant Term of Q(x)*P(x)^k #SEQ(binomial(2*k,k),k,p,100,2): SEQ:=proc(P,Q,x,N,r) local L,p,gu,mu,i: L:=CTtoSeq(P,Q,x,r*N+1): p:=5: gu:=[]: while p<=N do mu:=mods(convert([op(1..r*p,L)],`+`),p): gu:=[op(gu),[p,mu]]; p:=nextprime(p): od: gu,{seq(gu[i][2],i=1..nops(gu))} : end: #PaSu1(P,Q,x,p,r): inputs Laurent polynomials P,Q, (with integer coefficients)in the variable x, and a symbol p (denoting a general prime), #and a positive integer r, outputs #a simplified expression for SEQ(P,Q,x,N,r) PaSu1:=proc(P,Q,x,p,r) local m,POL,POL1,i: m:=-ldegree(P,x): if coeff(P,x,-m)<>1 then print(`The coefficient of the lowest power of`, P, `should have been 1 `): RETURN(FAIL): fi: POL:=expand(x^m*P): POL1:=add(coeff(POL^r,x,i)*x^i, i=0..r*m-1): POL1:=subs(x=x^p,POL1): simplify(x^m*Q*POL1/(POL-x^m)/x^(m*r*p)): end: #EvSumSepG(P,Q,x,y,r): Like PaSu1(P,Q,x,p,r) but with the rational function part seperated. Outputs it in the form #where y stands for 1/x^p: #[Expression with x^p, RationalFunction of x]. EvSumSepG:=proc(P,Q,x,y,r) local m,POL,POL1,i: option remember: m:=-ldegree(P,x): if coeff(P,x,-m)<>1 then print(`The coefficient of the lowest power of`, P, `should have been 1 `): RETURN(FAIL): fi: POL:=expand(x^m*P): POL1:=add(coeff(POL^r,x,i)*x^i, i=0..r*m-1): POL1:=subs(x=1/y,POL1): #Q*POL1/(POL-x^m)/x^(m*(r*p-1)): [expand(POL1*y^(m*r)), factor(x^m*Q/(POL-x^m))] end: #PaSu(P,Q,x,p,C,r): Like PaSu1(P,Q,x,p,r) but also inputs a symbol, C, where C(i) denotes the coefficient of the outputted rational function (the second part of the output) #the first part in a linear compination of terms of the form C(i*p-j, for specific integers i and j. Try: #PaSu((1/x+2+x),1,x,p,C,2); PaSu:=proc(P,Q,x,p,C,r) local y,gu,mu,R,gu1,i,j: gu:=EvSumSepG(P,Q,x,y,r): if gu=FAIL then RETURN(FAIL): fi: R:=gu[2]: gu:=gu[1]: mu:=0: for i from 0 to degree(gu,y) do gu1:=coeff(gu,y,i): for j from ldegree(gu1,x) to degree(gu1,x) do mu:=mu+C(i*p-j)*coeff(gu1,x,j): od: od: [mu,R]: end: #TheoG(P,Q,x,p,C,r): inputs Laurent polynomials P and Q, in the variable x, and symbols C and p, and a positive integer r,outputs a theorem #about the partial sum (up to p-1) mod p of the combinatorial sequence defined by C(n):=ConstantTerm (Q*P^n). #Try: #TheoG(1/x+2+x,1,x,p,C,r): TheoG:=proc(P,Q,x,p,C,r) local gu,A,i: gu:=PaSu(P,Q,x,p,C,r): if gu=FAIL then RETURN(FAIL): fi: print(``): print(`--------------------------------------------------------`): print(``): print(`Theorem: Define the following sequence: A(i) is the Constant term of the Laurent polynomial `): print(Q*P^i): print(`and for any prime p, let`): print(B(p)=Sum(A(i),i=0..r*p-1)): print(`Let the C-finite sequence C(n) be defined in terms of the generating function`): print(``): print(Sum(C(i)*x^i,i=0..infinity)=gu[2]): print(`then, modulo p, `): print(B(p)=gu[1]): print(``): end: #PaSu1check(P,Q,x,N,r): checks the output for PaSu1(P,Q,x,p,r) (or, equivalently, PaSu(P,Q,x,p,C,r), up to all primes<=N. #Try: #PaSu1check(1/x+2+x,1,x,100,2); PaSu1check:=proc(P,Q,x,N,r) local y,gu,lu,mu,i,j,mu1,p,coe1: gu:=EvSumSepG(P,Q,x,y,r): lu:=taylor(gu[2],x=0,N*r*degree(gu[1],y)+1): lu:=[seq(coeff(lu,x,i),i=1..N*r*degree(gu[1],y))]: mu:=[]: p:=5: while p<=N do mu1:=0: for i from 1 to degree(gu[1],y) do for j from 0 to degree(coeff(gu[1],y,i),x) do coe1:=coeff(coeff(gu[1],y,i),x,j): mu1:=mu1+coe1*lu[i*p-j]: od: od: mu:=[op(mu),[p, mods(mu1,p)]]: p:=nextprime(p): od: evalb(mu=SEQ(P,Q,x,N,r)[1]): end: #EvQPgOld(P,Q,x,p,C,r,d,n): Like PaSu(P,Q,x,p,C,r), but expresses C(i) as A quasi-polynomial of degree d, in n, rather than a sequence of coefficients of a rational function. #Try: #EvQPgOld((1/x+2+x),1,p,C,2,0,n); EvQPgOld:=proc(P,Q,x,p,C,r,d,n) local gu,R,mu: gu:=PaSu(P,Q,x,p,C,r): if gu=FAIL then RETURN(FAIL): fi: R:=gu[2]: mu:=RatToQP(R,x,300,d,n): if mu=FAIL then FAIL: else [gu[1],mu]: fi: end: #TheoQP(P,Q,x,p,r,d): inputs Laurent polynomials P and Q, in the variable x, and symbols C and p, and a positive integer r,outputs a theorem #about the partial sum (up to p-1) mod p of the combinatorial sequence defined by C(n):=ConstantTerm (Q*P^n). #expressed in terms of quasipolynomials oof n f degree<=d #Try: #TheoQP(1/x+2+x,1,x,p,2,0): TheoQP:=proc(P,Q,x,p,r,d) local gu,A,i,PER: gu:=CHZ(P,Q,x,r,d,p): if gu=FAIL then RETURN(FAIL): fi: PER:=nops(gu): print(``): print(`--------------------------------------------------------`): print(``): print(`Theorem: Define the following sequence: A(i) is the Constant term of the Laurent polynomial `): print(Q*P^i): print(`and for any prime p, let`): print(B(p)=Sum(A(i),i=0..r*p-1)): print(`then B(p) mod p equals to the following according to its remainder when divided by`, PER): for i from 1 to PER-1 do if igcd(i,PER)=1 then print(` if p mod `, PER, `equals `, i , `then, B(p) mod p equals `, gu[i], `mod p `): fi: od: end: #CHZ(P,Q,x,r,d,p): Like PaSu(P,Q,x,p,C,r), but finds the explicit quasi-polynomial ifit exists, representing the output #Try: #CHZ((1/x+2+x),1,x,3,0,p); CHZ:=proc(P,Q,x,r,d,p) local n,gu,R,mu,T1,T2,PER,y,i,j: option remember: gu:=EvSumSepG(P,Q,x,y,r): R:=gu[2]: gu:=gu[1]: mu:=RatToQP(R,x,300,d,n): if mu=FAIL then RETURN(FAIL): fi: PER:=nops(mu): T1[0]:=mu[PER]: for i from 1 to PER-1 do T1[i]:=mu[i]: od: for i from 0 to PER-1 do T2[i]:=add(coeff(gu,y,j)*subs(n=j*p, T1[j*i mod PER]),j=0..degree(gu,y)): od: [seq(T2[i],i=1..PER-1),T2[0]]: end: #CHZcheck(P,Q,x,r,d,N): Nunerically checks CHZ(P,Q,x,r,d,p) for all primes <=N. #Try: #CHZcheck((1/x+2+x),1,x,3,0,200); CHZcheck:=proc(P,Q,x,r,d,N) local p,gu,mu,lu,p0: gu:=CHZ(P,Q,x,r,d,p): if gu=FAIL then print(`There is probably no quasi-polynomial of degree`,d, ` representation, so there is nothing to check, use procedure PaSu1 to find the representation in terms of C-finite sequences `): print(`or else make d bigger. `): RETURN(FAIL): fi: mu:=[]: p0:=5: while p0<=N do mu:=[op(mu), [p0, mods(EvalQP(gu,p,p0),p0) ]]: p0:=nextprime(p0): od: lu:=SEQ(P,Q,x,N,r)[1]: evalb(mu=lu): end: