#GuessPol1(L,n,d): inputs a list of pairs [[input,output], ...] a variable name n
#and a pos. integer d, outputs a polynomial of degree d f such that f(input)=output
GuessPol1:=proc(L,n,d) local a,i,f,var,eq:

if nops(L)<d+1 then
 print(`underdetermined make`, d, `smaller . `):
  RETURN(FAIL):
fi:
#f:=a[0]+a[1]*n+...+a[d]*n^d:
var:={ seq(a[i],i=0..d)}:
f:=add(a[i]*n^i,i=0..d):
# f(L[i][1])=L[i][2]

eq:={seq( subs(n=L[i][1], f)=L[i][2], i=1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
fi:

subs(var,f):

end:

#GuessPol(L,n): inputs a list of pairs [input,output] and tries to guess
#a polynomial of degree <=nops(L)-3
GuessPol:=proc(L,n) local f,d:

for d from 0 to nops(L)-3 do
 f:=GuessPol1([op(1..d+3,L)] ,n,d):
  if f<>FAIL then
   f:=GuessPol1(L,n,d):
     if f<>FAIL then
       RETURN(f):
     fi:
   fi:
od:

print(`Too bad, it did not work out, generate more date, or bad news, it is not a polynomial`):
FAIL:

end:

 #Comp(n,k): inputs non-neg. integers n and k and outputs the SET of
 # compositions of n into k parts (vectors of non-neg. integers that add-up to n)
 #of length k
 Comp:=proc(n,k) local S,a1,S1, s1:
 option remember:
 
 if k<0 or n<0 then
   RETURN({}):
 fi:


  if n=0 then
    RETURN({[0$k]}):
  fi:

   if k=0 then
     if n=0 then
       RETURN({[]}):
     else
       RETURN({}):
     fi:
   fi:

 S:={}:

  for a1 from 0 to n do
     S1:=Comp(n-a1,k-1):
     S:=S union {seq([a1,op(s1)], s1 in S1)}:
   od:
S:

end:
 








    





















Help:=proc(): print(`Trunc(f,var,N), MRec(k,m,N), Coeff1(f,var,L), Mrec(k,m,N), Msq(k,N), M34(N) , Num(S,k,c) , Mab(a,b,n,K)`): end:


Num:=proc(S,k,c) local s,c1,su:
option remember:
if k=1 then
 if member(c,S) then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:

su:=0:

for s in S do
 c1:=c-s:
 su:=su+Num(S,k-1,c1):
od:

su:

end:



#Mab(a,b,n,K): the number of b by a magic rectangles with row sums a1*n and column sums b1*n as as a polynomial
#where a1=a/gcd(a,b) and b1=b/gcd(a,b)
#using K+1 data points, or the sequence if it can't find a polynomial
Mab:=proc(a,b,n,K) local gu,mu,i,a1,b1:

a1:=a/gcd(a,b):
b1:=b/gcd(a,b):

gu:= [seq([i,Num(Comp(a1*i,a),b,[(b1*i)$a])],i=0..K)]:

mu:=GuessPol(gu,n):

if mu<>FAIL then
 RETURN(mu):
else
RETURN([seq(Num(Comp(a1*i,a),b,[(b1*i)$a]),i=0..K)]):
fi:

end: