#John Kim
#Use whatever you like

Help:=proc(): print(` GP(L,x),GQP1(L,x,m), GQP(L,x), GFtoSeq(f,q,m), pmn(m,n), F(m,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:


#GFtoSeq(f,q,m): inputs an expression f in the variable q and 
#outputs the list of numbers L of length m, such that L[i] is 
#the coefficient of q^(i-1) in the Maclaurin series.

GFtoSeq:=proc(f,q,m) local i:
[seq(coeff(taylor(f,q=0,m),q,i),i=0..m-1)]:
end:


#pmn(m,n): inputs a pos. integer m, and outputs a quasipolynomial 
#of degree m and order lcm(1,..,m) that describes p_m(n).

pmn:=proc(m,n) local f,q,i:

f:=1/product(1-q^i,i=1..m):
GQP(GFtoSeq(f,q,(m+2)*lcm(seq(i,i=1..m))),n):

end:


#F(m,x): expresses trunc(x/m) as a quasi-polynomial in x of order m for integer x.

F:=proc(m,x) local L:

L:=[seq(trunc(i/m),i=1..6*m)]:
GQP(L,x):

end:


#T(n) = p_3(n)-trunc(n/4)*trunc(n+2/4)
#T(n) = [7/12+(1/3)*n+(1/12)*n^2, (1/12)*n^2+(1/3)*n, (1/12)*n^2+(1/3)*n+1/4, 
#        (1/12)*n^2+(1/3)*n+1/3, (1/12)*n^2+(1/3)*n+1/4, (1/12)*n^2+(1/3)*n]
#        -[(1/4)*n-1/4, (1/4)*n-1/2, (1/4)*n-3/4, (1/4)*n]
#        *[(1/4)*n-1/4, (1/4)*n+1/2, (1/4)*n+1/4, (1/4)*n]

#T(n) = T := [(1/48)*n^2+(11/24)*n+25/48, (1/48)*n^2+1/4+(1/3)*n, (1/48)*n^2+(11/24)*n+7/16, 
              (1/48)*n^2+(1/3)*n+1/3, (1/48)*n^2+(11/24)*n+3/16, (1/48)*n^2+1/4+(1/3)*n, 
              (1/48)*n^2+(11/24)*n+37/48, (1/48)*n^2+(1/3)*n, (1/48)*n^2+(11/24)*n+3/16, 
              (1/48)*n^2+7/12+(1/3)*n, (1/48)*n^2+(11/24)*n+7/16, (1/48)*n^2+(1/3)*n]