###################################################################### ##AFS.txt: Save this file as AFS.txt # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read `AFS.txt` # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Nov. 16, 2015 print(`Created: Nov. 16, 2015`): print(` This is AFS.txt `): print(`It the package that accompanies the article `): print(` Computerizing the Andrews-Fraenkel-Sellers Proofs on the Number of m-ary partitions mod m`): print(`(and doing MUCH more!) `): print(`by Shalosh B. Ekhad and Doron Zeilberger`): print(`and also available from Zeilberger's website`): print(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg/ .`): 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 Supporting procedures type ezra1();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`---------------------------------------`): print(`---------------------------------------`): print(`For a list of the Generalized procedures type ezraG();, 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(`---------------------------------------`): with(combinat): ezraS:=proc() if args=NULL then print(` The STORY procedures are: CapsuleGQv, CapsuleV `): print(``): else ezra(args): fi: end: ezraG:=proc() if args=NULL then print(` The Generalized procedures are: CapsuleG, EvalCapG, EvalCapGQ, CapsuleGQ, FEmG, ModMsG `): print(``): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(` The supporting procedures are: Atslan, AFSv, Bd, CheckCapsule, CheckCapsuleG, CheckCapsuleGQ, DZv, mSectR, mSectRg, RatToQP `): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: Capsule, CheckT1, EvalCap, FEm, ModMs `): print(` `): elif nops([args])=1 and op(1,[args])=AFSv then print(`AFSv(n,m): The original version of Theorem 2.1. of the Andrews-Fraenkel-Sellers Theorem 2.1 in Discrete Math 339(2016), 283-287.`): print(`Try:`): print(` AFSv(101,4); `): elif nops([args])=1 and op(1,[args])=Atslan then print(`Atslan(S,R,Sj,Rj,Sjm,Rjm): derives (once and for all) the Functional Equation satisfied by the j(mod m) part of a formal power series`): print( ` C(q), satisfying the Functional Equation `): print(` C(q)=S(q)+R(q)*C(q^m), where Sj, Rj are the respective j (mod m) parts. The output is the pair [A,B] such that `): print(` A and B are such that. `): print(` Cj(q)=A(q)+B(q)*Cj(q^m). Rjm is Rj(q^m) and Sjm is Rj(q^m). Try: `): print(` Atslan(S,R,Sj,Rj,Sjm,Rjm); `): elif nops([args])=1 and op(1,[args])=Bd then print(`Bd(n,d): inputs an integer n, and a base d, returns the base d list. Try:`): print(`Bd(10^10,7);`): elif nops([args])=1 and op(1,[args])=Capsule then print(`Capsule(R,q,m,i): Inputs a rational function R in the variable q, a positive integer m>1 and an integer `): print(`i, 0<=i1 and an integer `): print(` i, 0<=iq satisfies the functional equation `): print(` g(q)= S(q) g(q^m). BEFORE taking it mod m. Try: `): print(`FEm(1/(1-q),q,3,0);`): elif nops([args])=1 and op(1,[args])=FEmG then print(`FEmG(FE,q,m,i): Inputs a functional equation FE=[R,S], where R and S are rational functions of`): print(` of q, an intger m larger than 1, and an inteter 0<=inops(L)-2 then RETURN(FAIL): 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,n); 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 ##from Cfinite #CtoR(S,t), the rational function, in t, whose coefficients #are the members of the C-finite sequence S. For example, try: #CtoR([[1,1],[1,1]],t); CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L: if not (type(S,list) and nops(S)=2 and type(S[1],list) and type(S[2],list) and nops(S[1])=nops(S[2]) and type(t, symbol) ) then print(`Bad input`): RETURN(FAIL): fi: D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])): N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])): L1:=expand(D1*N1): L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1): f:=L1/D1: L:=degree(D1,t)+10: f1:=taylor(f,t=0,L+1): if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)): RETURN(FAIL): else RETURN(f): fi: end: #RtoC(R,t): Given a rational function R(t) finds the #C-finite description of its sequence of coefficients #For example, try: #RtoC(1/(1-t-t^2),t); RtoC:=proc(R,t) local S1,S2,D1,R1,d,f1,i,a: R1:=normal(R): D1:=denom(R1): d:=degree(D1,t): if degree(numer(R1),t)>=d then print(`The degree of the top must be less than the degree of the bottom`): RETURN(FAIL): fi: a:=coeff(D1,t,0): S2:=[seq(-coeff(D1,t,i)/a,i=1..degree(D1,t))]: f1:=taylor(R,t=0,d+1): S1:=[seq(coeff(f1,t,i),i=0..d-1)]: [S1,S2]: end: #GuessRec1(L,d): inputs a sequence L and tries to guess #a recurrence operator with constant cofficients of order d #satisfying it. It returns the initial d values and the operator # as a list. For example try: #GuessRec1([1,1,1,1,1,1],1); GuessRec1:=proc(L,d) local eq,var,a,i,n: if nops(L)<=2*d+2 then print(`The list must be of size >=`, 2*d+3 ): RETURN(FAIL): fi: var:={seq(a[i],i=1..d)}: eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): else RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]): fi: end: #GuessRec(L): inputs a sequence L and tries to guess #a recurrence operator with constant cofficients #satisfying it. It returns the initial values and the operator # as a list. For example try: #GuessRec([1,1,1,1,1,1]); GuessRec:=proc(L) local gu,d: for d from 1 to trunc(nops(L)/2)-2 do gu:=GuessRec1(L,d): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: #SeqFromRec(S,N): Inputs S=[INI,ope] #where INI is the list of initial conditions, a ope a list of #size L, say, and a recurrence operator ope, codes a list of #size L, finds the first N0 terms of the sequence satisfying #the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L). #For example, for the first 20 Fibonacci numbers, #try: SeqFromRec([[1,1],[1,1]],20); SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope: INI:=S[1]:ope:=S[2]: if not type(INI,list) or not type(ope,list) then print(`The first two arguments must be lists `): RETURN(FAIL): fi: L:=nops(INI): if nops(ope)<>L then print(`The first two arguments must be lists of the same size`): RETURN(FAIL): fi: if not type(N,integer) then print(`The third argument must be an integer`, L): RETURN(FAIL): fi: if N1 then print(`The constant term of R must be 1`): RETURN(FAIL): fi: lu:=taylor(R,q=0,N+1): lu:=add(coeff(lu,q,i)*q^i,i=0..N) : gu:=1: while ldegree(lu-1,q)<=N do gu:=gu*lu: gu:=add(coeff(gu,q,i)*q^i,i=0..N): lu:=subs(q=q^m,lu): od: [seq(coeff(gu,q,i),i=1..N)]: end: #FEm(R,q,m,i): Inputs a rational function R of q, an intger m larger than 1, and an inteter 0<=iq satisfies the functional equation #g(q)= S(q) g(q^m). Try: #FEm(1/(1-q),q,3,0); FEm:=proc(R,q,m,i) local Ri,lu: option remember: Ri:=mSectR(R,q,m,i): lu:=normal(Ri/subs(q=q^m,Ri)*R): lu:=convert(lu,parfrac): end: #Bd(n,d): inputs an integer n, and a base d, returns the base d list. Try: #Bd(10^10,7); Bd:=proc(n,d) local n1,m,gu: if n=0 then RETURN([]): elif n0 and m>1 and type(Ca,list) and nops(Ca)=m) then print(`Bad input`): RETURN(FAIL): fi: lu:=Bd(n1,m): mul(Ca[lu[i]+1],i=1..nops(lu)) mod m: end: #Capsule(R,q,m,i): Inputs a rational function R in the variable q, a positive integer m>1 and an integer #i, 0<=i0 then RETURN(FAIL): fi: a:=coeff(taylor(R,q=0,i+1),q,i) mod m: [a,[seq(coeff(lu,q,i1),i1=0..max(m-1,degree(lu,q) ) )] mod m]: end: #EvalCap(n1,m,Ca): inputs a positive integer n1, a modulo m, and a capsule Ca=[a,Capsule] #it outputs the value implied by the fast (log time) recurrence given by the capsule # Try: #EvalCap(100,3,[1,[1,2,0]]); EvalCap:=proc(n1,m,Ca) local a,C,INI,r1: option remember: INI:=Ca[1]: C:=Ca[2]: if n1<0 then RETURN(0): elif n1=0 then RETURN(INI): elif n10 then gu1:=[seq(EvalCap(n1,m,C),n1=0..N)]: gu2:=ModMs(R,q,m,(N+1)*m) mod m: gu2:=[seq(gu2[m*n1+i],n1=0..N)]: else gu1:=[seq(EvalCap(n1,m,C),n1=1..N)]: gu2:=ModMs(R,q,m,(N+1)*m) mod m: gu2:=[seq(gu2[m*n1+i],n1=1..N)]: fi: evalb(gu1=gu2): end: #CheckT1(m): Checks Theorem 1 in the paper CheckT1:=proc(m) local Ca, INI, C,i,x: if m mod 4=2 then print(`There is no theorem`): RETURN(FAIL): fi: if m mod 2 =1 then INI:=(m+1)/2: C:=[seq( op([binomial(i,2) mod m,binomial(i,2) mod m]),i=2..m-1)]: Ca:=[INI,C]: elif m mod 4 =0 then INI:=m/2: C:=[seq(trunc((i+2)*(i+4)*(2*i+3)/24) mod m ,i=0..m-3)]: C:=[op(C),seq(C[nops(C)-i+1],i=1..nops(C))]: Ca:=[INI,C]: fi: evalb(Capsule(1/(1-x)/(1-x^2),x,m,m-1)=Ca): end: #####################Start Generalized Portion for Functional Equations of the form f(q)=S(q)+R(q)*f(q^m) #ModMsG(S,R,q,m,N): Inputs rational functions S(q), R(q) of q , and a large positive integer N #finds the first N terms of the unique formal power series f(q) that satisfies the Functional Equation #f(q)=S(q)+R(q)*f(q^m). For example, try: #ModMsG(1,q/(1-q),q,3,100); ModMsG:=proc(S,R,q,m,N) local gu1,gu2, gu3, lu,i,mu: if S=0 then if subs(q=0,R)<>1 then print(`The constant term of R must be 1`): RETURN(FAIL): fi: fi: lu:=taylor(R,q=0,N+1): lu:=add(coeff(lu,q,i)*q^i,i=0..N) : mu:=taylor(S,q=0,N+1): mu:=add(coeff(mu,q,i)*q^i,i=0..N) : if S<>0 then gu1:=mu: else gu1:=subs(q=0,R): fi: gu2:=mu+lu*subs(q=q^m,gu1): gu2:=add(coeff(gu2,q,i)*q^i,i=0..N) : while gu1<>gu2 do gu3:=mu+lu*subs(q=q^m,gu2): gu3:=add(coeff(gu3,q,i)*q^i,i=0..N) : gu1:=gu2: gu2:=gu3: od: [seq(coeff(gu2,q,i),i=1..N)]: end: #mSectRg(R,t,m,i): Given a rational function R of t finds the #rational function, let's call it S(t) such that #the part of R with powers of the form t^(m*j+i) equals t^i*S(t^m). #Try: #mSectRg(t^2/(1-t-t^2),t,2,0); mSectRg:=proc(R,t,m,i) local R0,mana,gu1,gu2,j: mana:=quo(numer(R),denom(R),t): R0:=normal(R-mana): gu1:=add(coeff(mana,t,i+m*j)*t^j,j=0..trunc((degree(mana,t)-i)/m)): gu2:=mSectR(R0,t,m,i): normal(gu1+gu2): end: #Atslan(S,R,Sj,Rj,m): derives (once and for all) the Functional Equation satisfied by the j(mod m) part of a formal power series #C(q), satisfying the Functional Equation #C(q)=S(q)+R(q)*C(q^m), where Sj, Rj are the respective j (mod m) parts. The output is the pair [A,B] such that #A and B are such that. #Cj(q)=A(q)+B(q)*Cj(q^m). Rjm is Rj(q^m) and Sjm is Rj(q^m). Try: #Atslan(S,R,Sj,Rj,Sjm,Rjm): Atslan:=proc(S,R,Sj,Rj,Sjm,Rjm) local gu,Cj,Cjm: gu:=numer((Cj-Sj)/Rj -S-R*(Cjm-Sjm)/Rjm): gu:=solve(gu,Cj): normal([coeff(gu,Cjm,0), coeff(gu,Cjm,1)]): end: #FEmG(FE,q,m,j): Inputs a functional equation FE=[R,S], where # R and S are rational functions of q, an intger m larger than 1, and an inteter 0<=j0 then print(`Bad input`): RETURN(FAIL): fi: S:=FE[1]: R:=FE[2]: Rj:=mSectRg(R,q,m,j): Sj:=mSectRg(S,q,m,j): Rjm:=subs(q=q^m, Rj): Sjm:=subs(q=q^m, Sj): normal([-(R*Rj*Sjm-Rjm*Sj-S*Rj*Rjm)/Rjm, R*Rj/Rjm]): end: #CapsuleG(FE,q,m,i): Inputs am inhomog. Functional Equation, FE=[S,R], #where S and R are rational functions of the variable q, a positive integer m>1 and an integer #i, 0<=i0 then print(`Bad input`): RETURN(FAIL): fi: lu:=FEmG(FE,q,m,i): if subs(q=0,denom(lu[1]))=0 or subs(q=0,denom(lu[2]))=0 then print(`The functional equation needs some trivial change of dependent variable`): print(`not yet implemented`): RETURN(FAIL): fi: lu1:=lu[1]: lu2:=lu[2]: lu2:=convert(lu2,parfrac) mod m: if diff(denom(normal(lu2)),q)<>0 then RETURN(FAIL): fi: C:=[seq(coeff(lu2,q,i1),i1=0..max(m-1,degree(lu2,q) ) )] mod m: gu:=ModMsG(op(FE),q,m,nops(C)*m+i+4): gu:=[seq(gu[i+i1*m],i1=0..nops(C)-1)] mod m: [lu1,C,gu]: end: #EvalCapG(n1,m,Ca,q): inputs a positive integer n1, a modulo m, and a GENERALIZED capsule Ca=[A,Capsule,INI] #it outputs the value implied by the fast (log time) recurrence given by that generalized capsule # Try: #EvalCapG(100,3,[1/(1-q),[1,2,0],[1,2,0]],q); EvalCapG:=proc(n1,m,Ca,q) local A,C,INI,r1,ka,a: option remember: A:=Ca[1]: C:=Ca[2]: INI:=Ca[3]: if n1<0 then RETURN(0): elif n11 and an integer #i, 0<=im then print(`If j is larger then`, nops(C1)-1, `then C[j]=0`): fi: print(``): print(`d(n) modulo`, m, `can be computed VERY fast (in logarithmic time!) using the following recurrence, taken modulo`, m): K:=trunc((nops(C1)-1)/m): print(d(m*n+a)=add(C[a+j*m]*d(n-j),j=0..K)): print(``): print(`subject to the initial condition`): print(``): print(d(0)=INI): print(``): print(`Proof: `): print(`Since we are interested in`, d(n)=c(m*n+i), `Let's consider the generating function of d(n), let's call it G(q)`): print(``): print(G(q)= Sum(d(n)*q^n,n=0..infinity)): print(``): print(`Let our original generating function, given in the title, be called f(q) `): print(``): print(f(q)=Product(subs(q=q^(m^j),R),j=0..infinity)): print(``): print(`Obviously, f(q) saisfies the Functional Equation`): print(``): print(f(q)=R*f(q^m)): print(``): print(`Let us`, m, `-sect both sides and take the powers that are congruent to`, i, `modulo `, m): print(``): Ri:=mSectR(R,q,m,i): print(`It is readily seen that the part belonging to the powers that are`, i, `modulo `, m, `in the formal power series expansion of`, R): print(`can be written as`, q^i*R0(q^m) , `where R0 is the rational function`): print(``): print(Ri): print(``): print(`Exctacting the`, i, `mudulu `, m, `powers from both sides, we get`): print(``): print(q^i*G(q^m)=q^i*R0(q^m)*f(q^m) ): print(``): print(` Replacing `, q^m, `by `, q, `we get `): print(``): print(G(q)=R0(q)*f(q) ): print(`Hence `): print(f(q)=G(q)/R0(q) ): print(`Going back to the defining functional equation of f(q), this translates to a Functional Equation for G(q) `): print(G(q)/Ri=R*G(q^m)/subs(q=q^m,Ri)): print(`In other words `): print(``): print(G(q)=Ri*R/subs(q=q^m,Ri)*G(q^m)): print(``): lu:=convert(Ri*R/subs(q=q^m,Ri),parfrac): print(`Doing the partial fraction decomposition we get`): print(Ri*R/subs(q=q^m,Ri)=lu): print(`and taking it modulo`, m, ` we get `): lu:=lu mod m: print(``): print(Ri*R/subs(q=q^m,Ri)=lu): print(` `): print(`Hence, modulo`, m, `we have the functional equation `): print(``): print(G(q)=lu*G(q^m)): print(``): print(`and it is readily seen that this is equivalent to the statement of the theorem. QED `): print(``): print(`Just for fun the googol-th through googol+100-th terms of d(n) modulo`, m, `is `): print(seq( EvalCap(10^100+i1,m,[INI,C1]),i1=10^100..10^100+100)): print(``): print(`---------------------------------------------------------------------------------------------`): print(``): print(`This took`, time()-t0, `seconds, `): print(``): print(`---------------------------------------------------------------------------------------------`): print(``): end: #CapsuleGQ(FE,q,m,i,n,L): Like CapsuleG(FE,q,m,i) (q.v.) but the first component is a quasi-polynomial in n #(that is the explicit expression for the coeff. of q^n in A). L is a parameter for guessing #Try #CapsuleGQ([1,q/(1-q)],q,3,1,n,100); CapsuleGQ:=proc(FE,q,m,i,n,L) local gu,A,d: option remember: gu:=CapsuleG(FE,q,m,i): if gu=FAIL then RETURN(FAIL): fi: A:=gu[1]: d:=degree(denom(A),q): A:=RatToQP(A,q,L,d,n): if A=FAIL then FAIL: else [A,gu[2],gu[3]]: fi: end: #EvalCapGQ(n1,m,Ca,q,n): inputs a positive integer n1, a modulo m, and a GENERALIZED QUASI-POLYNOMIAL capsule Ca=[A,Capsule,INI] #it outputs the value implied by the fast (log time) recurrence given by that generalized capsule # Try: #EvalCapGQ(100,3,[[1],[1,2,0],[1,2,0]],q,n); EvalCapGQ:=proc(n1,m,Ca,q,n) local A,C,INI,r1,ka,a: option remember: A:=Ca[1]: C:=Ca[2]: INI:=Ca[3]: if n1<0 then RETURN(0): elif n1m then print(`If j is larger then`, nops(C1)-1, `then C[j]=0`): fi: if nops(INH)=1 then print(`Also let Q(n) be the polynomial`, INH[1]): else print(`Also let Q(n) be the quasi-polynomial of quasi-period`, nops(INH), ` such that`): print(`if n modulo`, nops(INH), ` is `, 0, ` then Q(n) equals `, INH[nops(INH)]): for i1 from 1 to nops(INH)-1 do print(`if n modulo`, nops(INH), ` is `, i1, ` then Q(n) equals `, INH[i1]): od: fi: print(``): print(`d(n) modulo`, m, `can be computed VERY fast (in logarithmic time!) using the following recurrence, taken modulo`, m): K:=trunc((nops(C1)-1)/m): print(d(m*n+a)=Q(m*n+a)+add(C[a+j*m]*d(n-j),j=0..K)): print(``): print(`subject to the initial condition`): print(seq(d(i1)=INI[i1+1],i1=0..nops(INI)-1)): print(`Just for fun the googol-th through googol+100-th terms of d(n) modulo`, m, `is `): print(seq( EvalCapGQ(10^100+i1,m,SHMOR,q,n),i1=10^100..10^100+100)): print(``): print(`---------------------------------------------------------------------------------------------`): print(``): print(`This took`, time()-t0, `seconds, `): print(``): print(`---------------------------------------------------------------------------------------------`): print(``): end: