print(`QuasiPol, created by Andrew Baxter. 2007.`): print(`Procedures Included: SeqDealer, GuessPol, GuessQuasiPol, QPEval`). print(`Type Help() to get more information in general. Use Help(ProcedureName) to get information on a specific procedure.`): print(`Given a sequence, GuessQuasiPol will attempt to guess a quasipolynomial function which generates the sequence.`): Help:=proc(): if nargs=0 then print(`A quasipolynomial is a piecewise-defined function F(n), where each value is determined by a polynomial in n and each condition is determined by the residue class of n modulo some integer m (called the condition modulus for f)`): print(`The main procedure in this package is GuessQuasiPol`): print(`Use Help(SeqDealer), Help(GuessPol), Help(GuessQuasiPol), Help(QPEval) for details about a specific procedure.`): elif args[1]=SeqDealer then print(`How to call it: SeqDealer(List,Modulus)`): print(`What it does: SeqDealer(L,M) distributes list L into A, a list of M lists.`): print(`Particulars: Element L[k*M+i] is assigned to A[i][k]. Assumes L begins with index 1.`): print(`Notes: The name comes from the idea of dealing cards from a deck L to M players seated around a table.`): elif args[1]=GuessPol then print(`How to call it: GuessPol(List,indeterminate)`): print(`What it does: GuessPol(L,x) returns either a polynomial in x which generates L or returns FAIL if it cannot find a polynomial of sufficiently small degree`): print(`Particulars: Base code comes from January 29, 2007 Experimental Math class. Maximum degree for returned polynomial is nops(L)-5`): print(`Notes: This is the engine which drives GuessQuasiPol. To improve the package, improve this procedure.`): elif args[1]=GuessQuasiPol then print(`How to call it: GuessQuasiPol(List,indeterminate, |list of guesses|)`): print(`What it does: GuessQuasiPol(L,x) returns either a quasipolynomial in x which generates list L or returns FAIL if no suitably small degree quasipolynomial can be found.`): print(`Particulars: Runs GuessQuasiPol1 for each possible condition modulus from 2 up until nops(L)/2.`): print(`Option: If a list (or set) of guesses for the condition modulus is included as a third argument, the procedure will try only these guesses`): print(`Notes: Using a list of guesses can accelerate the procedure, since there are fewer values to check`): elif args[1]=QPEval then print(`How to call it: QPEval(Quasipolynomial, Indeterminate, Argument)`): print(`What it does: Given quasipolynomial Q(x), QPEval(Q,x,N) returns Q(N)`): print(`Notes: The user can use QPEval in conjunction with the seq command to generate quasipolynomial sequences.`): fi: end: #SeqDealer(List,Modulus) distributes List into an array A of m=Modulus lists, where element List[mk+i] gets "dealt" (like cards) into A[i][k]. Assumes List is indexed starting with 1. SeqDealer:=proc(L::list,m::posint) local A, i, n: for i from 0 to m-1 do A[i]:=[]: od: for n from 1 to nops(L) do i:=modp(n,m): A[i]:= [op(A[i]),L[n]]: od: A[m]:=A[0]: [A[j]$j=1..m]: end: #GuessPol(L,x): inputs a sequence of numbers #and outputs a polynomial P of x, of degree nops(L)-5 #such that L[i]=P(i), for i=1..nops(L) GuessPol:=proc(L,x) local d,P,a,var,eq: d:=nops(L)-5: P:=add(a[i]*x^i,i=0..d): var:={seq(a[i],i=0..d)}: eq:={seq(subs(x=i,P)=L[i],i=1..nops(L))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: subs(var,P) end: #GuessQuasiPol(L,x) tries to guess a quasi-polynomial (in x) for the terms of L, including guessing the condition modulus. GuessQuasiPol:=proc(L,x) local L1,P,m,r,flag,CMList,i: CMList:=[i$i=2..nops(L)/2]: if nargs=3 then CMList:=args[3]: fi: for m in CMList do L1:=SeqDealer(L,m): flag:=GuessPol(L1[1],x): for r from 2 to m while flag<>FAIL do P[r-1]:=subs(x=1+(x-(r-1))/m,flag): flag:=GuessPol(L1[r],x): od: if flag<>FAIL then return [seq(P[i],i=1..m-1),subs(x=x/m, flag)]: fi: od: FAIL: end: QPEval:=proc(Q,x,N) local m,r, k: m:=nops(Q): r:=modp(N,m): k:=(N-r)/m: if r=0 then r:=m: fi: subs(x=N,Q[r]): end: #Generating Function Tools #These may be duplicated elsewhere in Maple, especially GenFuncExp. #GenFuncExp(GF,q,N) given generating function GF(q), returns the first N coefficients, starting with q^1. GenFuncExp:=proc(GF,q,N) local P: P:=convert(taylor(GF,q,N+1),polynom): [seq(coeff(P,q,i),i=1..N)]: end: #QPGF(A,B,q) returns 1/((1-q^A[1])^B[1] * (1-q^A[2])^B[2] ... (1-q^A[k])^B[k]) #Generating Functions of this form are known to have quasipolynomial generating functions. QPGF:=proc(A::list,B::list,q) local P: if nops(A)<>nops(B) then ERROR(`First two arguments must have same length`): fi: P:=mul((1-q^A[i])^B[i],i=1..nops(A)): 1/P: end: