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: