Help:=proc(): print(` GP(L,x),GQP1(L,x,m), GQP(L,x) `): end:
Di:=proc(L) local i:
[ seq( L[i+1]-L[i],i=1..nops(L)-1)]:
end:
#GP(L,x): The unique polynomial P(x)
#expressed in terms of binomial(x,k)
#such that P(i)=L[i+1] for i from 0 to nops(L)-1
GP:=proc(L,x) local i,L1,P:
L1:=L:
P:=0:
for i from 1 to nops(L)-2 do
P:=P+L1[1]*binomial(x,i-1):
L1:=Di(L1):
if convert(L1,set)={0} then
RETURN(factor(expand( subs(x=x-1,P)))):
fi:
od:
FAIL:
end:
#GQP1(L,x,m): inputs a list of numbers L, a variable x
#and A pos. integer m, and outputs a quasi-poly of order m,
#let's call it Q
#fitting L such that if n is i mod m then L[n] equals
#Q[i](n)
GQP1:=proc(L,x,m) local Q,i,j,L1,P1:
Q:=[]:
for i from 1 to m do
L1:=[ seq(L[i+m*j],j=0..trunc((nops(L)-i)/m) ) ]:
P1:=GP(L1,x):
if P1=FAIL then
RETURN(FAIL):
else
P1:=expand(subs(x=(x-i)/m+1,P1)):
Q:=[ op(Q), P1]:
fi:
od:
Q:
end:
#GQP(L,x): inputs a list of numbers and tries to guess
# a quasi-polynomial fitting the data, or returns FAIL
GQP:=proc(L,x) local m, Ben:
for m from 1 to nops(L)/2 do
Ben:=GQP1(L,x,m):
if Ben<>FAIL then
RETURN(Ben):
fi:
od:
FAIL:
end: