#C6.txt: Feb. 11, 2019

Help:=proc(): print(` NuM(A,k,c) , Delt(L) , GFpower(x,k) , GFpower(x,k) `): end:

 #The generating function of the discrete function n^k
 GFpower:=proc(x,k)
 option remember:
 if k=0 then
     RETURN(1/(1-x)):
 else
 normal(x*diff(GFpower(x,k-1),x)):
 fi:
 end:

 #Delt(L): inputs a list L and outputs a list of L[i+1]-L[i]
 Delt:=proc(L) local i: [seq(L[i+1]-L[i],i=1..nops(L)-1)]: end:


 #NuM(A,k,c):The number of vectors of length k whose components are
 #members of A that add-up to c. For example the number of 3 by 3 magic squares
 #with row sums and column sums all equal to 5 is NuM(Comp(5,3),3,[5,5,5]);
 NuM:=proc(A,k,c) local r,a,su:
 option remember:
 if not (type(c,list) and type (A,set)) then
  print(`Bad input`):
   RETURN(FAIL):
 fi:

 r:=nops(c):

 if k=0 then
  if c=[0$r] then
   RETURN(1):
   else
    RETURN(0):
   fi:
 fi:
 su:=0:

  for a in A do
   su:=su+NuM(A,k-1,c-a):
  od:
 
 su:
 end:



##stuff from C4.txt and C5.txt
#C5.txt: Continuing the polynomial anstaz
Help5:=proc(): print(` PreMag(a,b,n), IsMagic(A,k) , MagRec(a,b,n)  `): end:
###stuff from C4.txt

Help4:=proc(): print(`GuessPol1(L,n,d) , GuessPol(L,n), Comp(n,k)  `): end:

#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:

#end of C4.txt

   #PreMag(a,b,n):inputs non-neg. integers n,a,b, and outputs the set of a by b pre-magic 
   #rectangle with row sums n
   PreMag:=proc(a,b,n) local  S,J,s,j:
   option remember:
   if a=0 then
     RETURN({[]}):
   fi:
    
 J:=Comp(n,b):

   S:=PreMag(a-1,b,n):

 {seq(seq([op(s),j], s in S), j in J)}:
  
  end:
 
 #IsMagic(A,k): inputs a matrix A and outputs true iff
 #all the columns of A add-up to k
 IsMagic:=proc(A,k) local i,j:
   evalb({seq(add(A[i][j],i=1..nops(A)),j=1..nops(A[1]) ) }= {k}):
 end:
   
  
 #MagRec(a,b,n): the set of a by b matrices with non-neg. entries such that
 #every row sum is n and every column sum a*n/b  
 MagRec:=proc(a,b,n) local A,G,c,a1:
   if not type(a*n/b,integer) then
     RETURN({}):
   fi:

   c:=a*n/b:
  
  A:=PreMag(a,b,n):
  
  G:={}:

  for a1 in A do
   if IsMagic(a1,c) then
    G:=G union {a1}:
   fi:
 od:

 G:
end:

#end C5.txt