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: