######################################################################
##MagicRectangles.txt: Save this file as  MagicRectangles.txt        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read MagicRectangles.txt<Enter>                    #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#and his Math640 (Fall 2019) class                                   #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Feb. 2019
 
print(`Created:  Feb. 2019`):
print(` This is MagicRectangles.txt `):
print(`It is the package that accompanies the article `):

print(` Explicit Polynomial Expressions for the Enumeration of 3 by 4 and 3 by 5 Semi-Magic Rectangles`):
print(` By Shalosh B. Ekhad and StudentsOfDrZ ExpMathClass`):

print(` Explicit Expressions for the Number of Semi-Magic Rectangles`):
print(`by StudentsOfDrZ ExpMathClass and Shalosh B. Ekhad`):
print(`published in the Personal Journal of Shalosh B. Ekhad and Doron Zeilberger`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Mab , SipurMR`):
 print(` `):



elif nops([args])=1 and op(1,[args])=Mab then
print(`Mab(a,b,n): the number of a by b magic rectangles with row sums b1*n and column sums a1*n as as a polynomial`):
print(`where a1=a/gcd(a,b) and b1=b/gcd(a,b). Try:`):
print(`Mab(3,3,n);`):

elif nops([args])=1 and op(1,[args])=SipurMR then
print(`#SipurMR(L,n,M,N): inputs a list of non-negative integers and outputs explicit expressions and the first M terms `):
print(`starting at n=1, as well the value at n=N`):
print(`the number of semi-magic rectangles of size axb for a from 1 to nops(L) and b from a to L[a]. Try:`):
print(`SipurMR([0,2,3,4],n,20,1000);`):

else
print(`There is no ezra for`,args):
fi:
 
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(factor(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)`): 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): the number of a by b magic rectangles with row sums b1*n and column sums a1*n as as a polynomial
#where a1=a/gcd(a,b) and b1=b/gcd(a,b). Try:
#Mab(3,3,n);
Mab:=proc(a,b,n) local gu,mu,i,a1,b1,K:

K:=(a-1)*(b-1)+2:
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:


##SipurMR(L,n,M,N): inputs a list of non-negative integers and outputs explicit expressions and the first M terms 
#starting at n=1, as well the value at n=N
#the number of semi-magic rectangles of size axb for a from 1 to nops(L) and b from a to L[a]. Try:
#SipurMR([0,2,3,4],n,20,1000);
SipurMR:=proc(L,n,M,N) local gu,a,b,co,n1,a1,b1:

co:=0:

print(`Explicit Polynomial Expressions for Certain Semi-Magic Rectanlges`):
print(``):
print(`By Shlosh B. Ekhad `):
print(``):
print(`In this article we will give polynomial expressions in n, for the number of semi-magic a by b rectangles`):
print(``):
for a from 1 to nops(L) do
 for b from a to L[a] do
  co:=co+1:
  gu:=Mab(a,b,n):
   a1:=a/gcd(a,b):
   b1:=b/gcd(a,b):
 print(``):
print(`----------------------------------------------`):
  print(``):
  print(`Theorem number`, co, `: The number of`, a, ` by ` , b, `rectangles with row-sums`, b1*n, `and column-sums`, a1*n, `equals the`):
  print(`the following polynomial in`, n, `of degree `, (a-1)*(b-1)):
  print(``):
   print(gu):
  print(``):
  print(`and in Maple format`):
   print(``):
  lprint(gu):
    print(``):
 print(`For the sake of the OEIS, the first`, M, `terms, starting at`, n, `= 1  are `):
    print(``):
  print([seq(subs(n=n1,gu),n1=1..M)]):
      print(``):
 print(`Just for fun, the number of`, a, ` by `, b, `semi-magic rectangles with row-sums`,a1*N,`and column-sums`, b1*N, ` equals`):
      print(``):
  lprint(subs(n=N,gu)):
       print(``):
 od:
od:

print(``):
print(`This ends this web-book that took`, time(), `seconds to generate. `):
print(``):
end: