######################################################################
##SALIKHOVpiHuman.txt: Save this file as  SALIKHOVpiHuman.txt        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  SALIKHOVpiHuman.txt<Enter>                   #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Dec.  2019- Jan. 2020

print(`Created:  Dec. 2019- Jan. 2020`):
print(` This is SALIKHOVpiHuman.txt `):
print(`It is  one of the Maple packages that accompanies the article `):
print(`   The Irrationality Measure of Pi is at most 7.103205334137... `):
print(`by  Doron Zeilberger and Wadim Zudiln`):
print(`available from the arxiv and from `):
print(` http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/pimeas.html `):
print(``):
print(`It only handles the human claims in Part II of the paper`):
print(``):
print(`For the Maple package accompanying part I see the Maple package SALIKHOVpi.txt`):
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 procedures are:  CheckL2G,  CheckL2pG, CheckL2v, Comps, CompsG, IsGoodpj, Mjn, OrdpTmn, Pn, rRangep, `):
 print(`TayR1modp, TayRfastP, Tmn `):
else
 ezra(args):
fi:
end:

ezra:=proc():

if args=NULL then
 print(`The main procedures are: Alist, AlistNegPre, AlistNeg, Anj,Anj5,  Anpj, `):
 print(`CheckEq5, CheckEq5a,CheckEq5b,  ChecEq5A, ChecEq5B,  CheckEq6A, CheckEq6B, CheckEq7, CheckEq7A, CheckEq7B,`):
 print(`CheckEq7C, CheckEq7D,  Checkj0, CheckL1,  CheckL2,`):
 print(` CheckL2A, CheckL2Apn, CheckL3, CheckL4, CheckL5, CheckL6, CheckP1, Ln, Np,  Pn, Pnx, Refusinks11, Refusinks1, RefusinksB `):
 print(`RefusinksZ1, RefusinksZ, Rnx, RnxS, TayR, TayRfast, WtEq5S , Znmj `):
 print(` `):


elif nops([args])=1 and op(1,[args])=Alist then
print(`Alist(n): the list of length 3*n+1 consisting of the A_j from j=0 to 3*n. Try:`):
print(`Alist(5);`):


elif nops([args])=1 and op(1,[args])=AlistNeg then
print(`AlistNeg(n): The list [A[-1], ..., A[-(4*n-1)] as given by Equation 6. Try:`):
print(`AlistNeg(3);`):

elif nops([args])=1 and op(1,[args])=AlistNegPre then
print(`AlistNegPre(n): The list of coefficients  of (x+5)^(k-1) in Pn(n,x) from k=1 to k=4*n-1. Try:`):
print(`AlistNegPre(3);`):


elif nops([args])=1 and op(1,[args])=Anj then
print(`Anj(n,j): the coefficient A_j, of 1/(x+5)^(j+1) in the partial fraction expansion of R(x).  Try:`):
print(`Anj(5,3);`):

elif nops([args])=1 and op(1,[args])=Anj5 then
print(`Anj5(n,j): The formula for A_j given in Eq(5)`):
print(`Try:`):
print(`Anj5(3,0);`):

elif nops([args])=1 and op(1,[args])=Anpj then
print(`Anpj(n,p,j): the A_j mod p in the n case`):
print(`Anpj(12,23,0);`):

elif nops([args])=1 and op(1,[args])=CheckEq2 then
print(`CheckEq2(N): Checks that the two expressions for R_n(x) given in Equation (2) are correct for n from 1 to N do`):
print(`CheckEq2(20);`):

elif nops([args])=1 and op(1,[args])=CheckEq5 then
print(`CheckEq5(n): checks equation (5) of the ZZ paper for n. Try:`):
print(`CheckEq5(4);`):


elif nops([args])=1 and op(1,[args])=CheckEq5a then
print(`CheckEq5a(n): checks the penultimate line of equation (5) of the ZZ paper for n. Try:`):
print(`CheckEq5a(4);`):


elif nops([args])=1 and op(1,[args])=CheckEq5b then
print(`CheckEq5b(n): checks the secondline of equation (5) of the ZZ paper for n. Try:`):
print(`CheckEq5b(4);`):


elif nops([args])=1 and op(1,[args])=CheckEq5A then
print(`CheckEq5A(n): checks the last line of p.5. Try:`):
print(`CheckEq5A(4);`):


elif nops([args])=1 and op(1,[args])=CheckEq5B then
print(`CheckEq5B(n): checks the first line of p.6. Try:`):
print(`CheckEq5B(4);`):


elif nops([args])=1 and op(1,[args])=CheckEq5C then
print(`CheckEq5C(n): checks the equation in  line 6 of p.6. Try:`):
print(`CheckEq5C(3);`):


elif nops([args])=1 and op(1,[args])=CheckEq6 then
print(`CheckEq6(n): checks the Eq. (6) of the ZZ paper for n. Try:`):
print(`CheckEq6(2);`):

elif nops([args])=1 and op(1,[args])=CheckEq6A then
print(`CheckEq6A(n): checks the correctness of the expression on line 4 for n and all j in {-4*n,..., 3*n}. Try:`):
print(`CheckEq6A(10);`):

elif nops([args])=1 and op(1,[args])=CheckEq6B then
print(`CheckEq6B(N,K): checks that the p-order of N! is trunc(N/p) if p>sqrt(N) up to K primes. Try:`):
print(`CheckEq6B(10,10); `):

elif nops([args])=1 and op(1,[args])=CheckEq7 then
print(`CheckEq7(N,K): Checks Equation (7) in the ZZ paper for all n<=N, n1<=2*n,  and the first K primes. Try:`):
print(`CheckEq7(10,10);`):

elif nops([args])=1 and op(1,[args])=CheckEq7A then
print(`CheckEq7A(N): Checks the alternative expressions for Z(n,m,j) given in the last two lines of p. 7 for`):
print(`n,m<=N and j between 0 and 3*n. Try:`):
print(`CheckEq7A(10);`):

elif nops([args])=1 and op(1,[args])=CheckEq7B then
print(`CheckEq7B(K): Checks  the identity in the third line of  p. 8 for 0<=M<=N<=K. Try`):
print(`CheckEq7B(10);`):


elif nops([args])=1 and op(1,[args])=CheckEq7C then
print(`CheckEq7C(K): Checks  that the sum the 6th line of  p. 8 for 0<=M,N<=K, 0<=j<=2*M. equals the coefficient of t^(2*N)`):
print(`in the polynomial in the 8th line Try`):
print(`CheckEq7C(10);`):

elif nops([args])=1 and op(1,[args])=CheckEq7D then
print(`CheckEq7D(N): Checks the statement about the p-adic order of (2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)!`):
print(`for all0<=n1<=n2<=2*n, p in P_n, and  j between -(4*n-1) to 3*n. Try:`):
print(`CheckEq7D(10);`):

elif nops([args])=1 and op(1,[args])=Checkj0 then
print(`Checkj0(N):checks Lemma2 for j=0 and all n<=N. Try:`):
print(`Checkj0(1000);`):

elif nops([args])=1 and op(1,[args])=CheckL1 then
print(`CheckL1(n): empirically checks Lemma 1 in the Zeilberger-Zudilin paper for n. Try: `):
print( `CheckL1(5); `):

elif nops([args])=1 and op(1,[args])=CheckL2 then
print(`CheckL2(n): checks Lemma 2. Try:`):
print(`CheckL2(10);`):

elif nops([args])=1 and op(1,[args])=CheckL2A then
print(`CheckL2A(n): checks the statement in p. 7, lines 4-7 of`):
print(` http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/pimeas.pdf for`):
print(`n, all primes in Pn(n) and 0<=n0,n1,n2<=2*n. Try:`):
print(`CheckL2A(10);`):

elif nops([args])=1 and op(1,[args])=CheckL2Apn then
print(`CheckL2Apn(p,n): Checks Lemma2A for a specific p and n Try:`):
print(`CheckL2Apn(7,9);`):

elif nops([args])=1 and op(1,[args])=CheckL2G then
print(`CheckL2G(N): checks Lemma 2 for all n from 1 to N do`):
print(`CheckL2G(100);`):

elif nops([args])=1 and op(1,[args])=CheckL2pG then
print(`CheckL2pG(p,G): Inputs a prime p, and a positive integer G, checks the generalized Lemma 2 for all n in NpG(p,G). Try:`):
print(`CheckL2pG(17,trunc(17^2/3));`):

elif nops([args])=1 and op(1,[args])=CheckL2v then
print(`CheckL2v(n): checks Lemma 2. Verbose version. Try:`):
print(`CheckL2v(10);`):


elif nops([args])=1 and op(1,[args])=CheckL3 then
print(`CheckL3(N): checks Lemma 3 on page 9 for n from 1 to N. Try:`):
print(`CheckL3(10);`):


elif nops([args])=1 and op(1,[args])=CheckL4 then
print(`CheckL4(N): checks Lemma 4 on page 10 for n from 1 to N. Try:`):
print(`CheckL4(20);`):


elif nops([args])=1 and op(1,[args])=CheckL5 then
print(`CheckL5(N): checks Lemma 5 on page 10 for n between 1 and N. Try:`):
print(`CheckL5(30);`):


elif nops([args])=1 and op(1,[args])=CheckL6 then
print(`CheckL6(N): checks Lemma 6 on page 11 for n between 1 and N. Try:`):
print(`CheckL6(10);`):


elif nops([args])=1 and op(1,[args])=CheckP1 then
print(`CheckP1(N): checks Proposition 1 on page 12 for n between 1 and N. Try:`):
print(`CheckP1(10);`):

elif nops([args])=1 and op(1,[args])=CheckOriL2 then
print(`CheckOriL2(n): checks the original Lemma 2. Try:`):
print(`CheckOriL2(10);`):

elif nops([args])=1 and op(1,[args])=Comps then
print(`Comps(N,A,r): The set of r-tuples of non-negative integers <=A that add-up to N. Try: `):
print(`Comps(5,3,3); `):

elif nops([args])=1 and op(1,[args])=CompsG then
print(`CompsG(N,A): Inputs a positive ineger N, a list  A, or length r, say ,of non-negative integers outputs`):
print(`all the vectors of length r [v1, ..., vr] such that v[i]<=A[i] and that add-up to N`):
print(`CompsG(5,[2,2,3]);`):

elif nops([args])=1 and op(1,[args])=IsGoodpj then
print(`IsGoodpj(p,j,n):Is p divisible by T(m) for all m in M_j?Try:`):
print(`IsGoodpjn(5,3,6);`):

elif nops([args])=1 and op(1,[args])=Ln then
print(`Ln(n):L_n defined in the statement of Lemma 3 on p. 9. Try:`):
print(`Ln(10);`):

elif nops([args])=1 and op(1,[args])=Mjn then
print(`Mjn(j,n): the set M_j defined on line 6 from the bottom of page 5. Try:`):
print(`Mjn(3,3);`):


elif nops([args])=1 and op(1,[args])=Np then
print(`Np(p): Inputs a prime p, larger than 5, outputs the set of n such that n<p^2/3 and 1/2<={n/p}<2/3*p. Try:`):
print(`Np(17);`):

elif nops([args])=1 and op(1,[args])=OrdpTmn then
print(`OrdpTmn(p,m,n): The lower bound for the p-adic order of Tmn(m,n). Try:`):
print(`OrdpTmn(11,[0,4,5,3,2,1],30);`):

elif nops([args])=1 and op(1,[args])=Pn then
print(`Pn(n): the set of primes defined in Lemma2. Try: `):
print(`Pn(11);`):

elif nops([args])=1 and op(1,[args])=Pnx then
print(`Pnx(n,x): The polynomial part of R_n(x). Try:`):
print(`Pnx(10,x);`):

elif nops([args])=1 and op(1,[args])=Refusinks11 then
print(`Refusinks11(n): inputs a positive integer n, and  outputs the set of pairs [p,Listp], where Liptp is the list of n1 where 0<=n1<=n`):
print(` such that binomial(2*n,n1) is NOT divisible by p.`):
print(`Try Refusinks11(10);`):



elif nops([args])=1 and op(1,[args])=Refusinks1 then
print(`Refusinks11(n): inputs a positive integer n, and  outputs the set of pairs [p,Listp], where Liptp is the list of n1 where 0<=n1<=2*n`):
print(` such that binomial(2*n,n1) is NOT divisible by p.`):
print(`Try Refusinks11(10);`):



elif nops([args])=1 and op(1,[args])=RefusinksB then
print(`RefusinksB(n): inputs a positive integer n, and  outputs the set of pairs [p,Listm], where Listm is the set of n1+n2 such that`):
print(`  binomial(2*n,n1)*binomial(2*n,n2) is NOT divisible by p.`):
print(`Try RefusinksB(10);`):



elif nops([args])=1 and op(1,[args])=RefusinksZ then
print(`RefuskinksZ(n): Inputs a positive integer n and output the list of pairs [p,Lisp] where Listp is the list`):
print(`of m such that Z(n,m,i*p) is not divisibile by p for at least one i. Try`):
print(`RefuskinsZ(20);`):

elif nops([args])=1 and op(1,[args])=RefusinksZ1 then
print(`RefuskinksZ1(n,i): Inputs a positive integer n and output the list of pairs [p,Lisp] where Listp is the list`):
print(`of m such that Z(n,m,i*p) is not divisibile by p. Try`):
print(`RefuskinsZ1(20,0);`):

elif nops([args])=1 and op(1,[args])=rRangep then
print(`rRangep(p): Inputs a prime p, outouts the set of integers between p/2 and 2*3/p. Try;`):
print(`rRangep(11);`):


elif nops([args])=1 and op(1,[args])=Rnx then
print(`Rnx(n,x): The rational function R_n(x) defined in the first part of Equation (2) of the ZZ paper. Try:`):
print(`Rnx(3,x);`):


elif nops([args])=1 and op(1,[args])=RnxS then
print(`RnxS(n,x): The rational function R_n(x) defined in the second part of Equation (2) (not using complex numbers) of the ZZ paper. Try:`):
print(`RnxS(3,x);`):

elif nops([args])=1 and op(1,[args])=TayR then
print(`TayR(n,u): the Taylor expansion of the rational function 5*(u-5)^(2*n)*(u^4-20*u^3+156*u^2-560*u+800)^(2*n)/(10-u)^(3*n+1) up to 3n+1 terms. Try:`):
print(`TayR(5,u);`):

elif nops([args])=1 and op(1,[args])=TayR1modp then
print(`TayR1modp(u,p,K): the Taylor expansion of `):
print(`R1=(u^4-20*u^3+156*u^2-560*u+800)^2*(u-5)^2/(10-u)^3 up to u^(K+1) mod p. Try:`):
print(`TayR1modp(u,11,30);`):

elif nops([args])=1 and op(1,[args])=TayRfast then
print(`TayRfast(n,u,p): Fast version of the Taylor expansion of`):
print(`5*(u-5)^(2*n)*(u^4-20*u^3+156*u^2-560*u+800)^(2*n)/(10-u)^(3*n+1) mod p`):
print(`up to 3*n+1 terms. Try:`):
print(`TayRfast(20,u,7);`):

elif nops([args])=1 and op(1,[args])=TayRfastP then
print(`TayRfastP(u,p,n,K): The first K coefficients of the Taylor exapnsion of`):
print(`R1^n where`):
print(`R1=(u^4-20*u^3+156*u^2-560*u+800)^2*(u-5)^2/(10-u)^3 up to u^(K+1) mod p. Try:`):
print(`TayRfastP(u,7,20,30);`):

elif nops([args])=1 and op(1,[args])=TayRmodp then
print(`TayRmodp(n,u,p): the Taylor expansion of the rational function 5*(u-5)^(2*n)*(u^4-20*u^3+156*u^2-560*u+800)^(2*n)/(10-u)^(3*n+1) up to 3n+1 terms`):
print(`modlulo p . It is given in two parts`):
print(`Try:`):
print(`TayRmodp(10,u,13);`):

elif nops([args])=1 and op(1,[args])=Tmn then
print(`Tmn(m,n): the quantity T(m) defined on line 3 from the bottom of page 5. m is a six-tuple of non-negative intgers and n is a pos. integer. Try: `):
print(`Tmn([0,1,2,3,4,5],5);`):


elif nops([args])=1 and op(1,[args])=WtEq5 then
print(`WtEq5(m,n): Given  a six-tuple m and  an integer  n,outputs the summand in the second line of equation `):
print(`Try:`):
print(`WtEq([0,0,0,0,0,0],4); `):

elif nops([args])=1 and op(1,[args])=WtEq5S then
print(`WtEq5S(m,n): Given  symbols m and n,outputs the symbolic expression for the summand of the second line of Eq. (5)`):
print(`Try:`):
print(`WtEq(m,n):`):

elif nops([args])=1 and op(1,[args])=Znmj then
print(`Znmj(n,m,j): The discrete function Z(n,m,j) defined in line 8 of p. 7, Try:`):
print(`Znmj(4,2,1);`):


elif nops([args])=1 and op(1,[args])=Z1nmj then
print(`Z1nmj(n,m,j): The alternative expression for Z(n,m,j) at the penultimate line of p. 7. Try:`):
print(`Z1nmj(4,2,1);`):


elif nops([args])=1 and op(1,[args])=Z2nmj then
print(`Z2nmj(n,m,j): The alternative exression for Z(n,m,j) at the last line of p. 7. Try:`):
print(`Z2nmj(4,2,1);`):


else
print(`There is no ezra for`,args):
fi:
 
end:



#AlistOld(n): the list of length 3*n+1 consisting of the A_j from j=0 to 3*n. Try:
#AlistOld(5);
AlistOld:=proc(n) local R,x,i,j,gu,gu1, A:

R:=5*x^(2*n)*(x^4+6*x^2+25)^(2*n)/(5+x)^(3*n+1)/(5-x)^(3*n+1):

gu:=convert(R,parfrac):


for i from 1 to nops(gu) do
 gu1:=op(i,gu):
 
 if degree(denom(gu1),x)>0 then
   j:=degree(denom(gu1),x)-1:
 
   A[j]:=numer(gu1):
 fi:
od:

[seq(A[j],j=0..3*n)]:

end:

#CheckL1(n): empirically checks Lemma 1 in the Zeilberger-Zudilin paper for n. Try
#CheckL1(5);
CheckL1:=proc(n) local R,x,i,j,gu,gu1, A:

R:=5*x^(2*n)*(x^4+6*x^2+25)^(2*n)/(5+x)^(3*n+1)/(5-x)^(3*n+1):

gu:=convert(R,parfrac):


for i from 1 to nops(gu) do
 gu1:=op(i,gu):
 
 if degree(denom(gu1),x)>0 then
   j:=degree(denom(gu1),x)-1:
 
   A[j]:=numer(gu1):
 fi:
od:

evalb({seq(type(2*A[j]/5^j/2^(trunc((5*n+3*j)/2)),integer),j=0..3*n)}={true}):

end:


#Comps(N,A,r): The set of r-tuples of non-negative integers <=A that add-up to N. Try
#Comps(5,3,3);
Comps:=proc(N,A,r) local gu,i,lu,lu1:
option remember:

if N<0 then
 RETURN({}):
fi:

if N=0 then
 RETURN({[0$r]}):
fi:
if r=1  then
  if N<=A then
   RETURN({[N]}):
  else
   RETURN({}):
  fi:
fi:

gu:={}:

for i from 0 to A do
 lu:=Comps(N-i,A,r-1):
 gu:=gu union {seq([op(lu1),i],lu1 in lu)}:
od:

gu:

end:


#Mjn(j,n): the set M_j defined on line 6 from the bottom of page 5. Try:
#Mjn(3,3);
Mjn:=proc(j,n)
option remember:
CompsG(3*n-j,[3*n-j,2*n$5]):
end:

#Tmn(m,n): the quantity T(m) defined on line 3 from the bottom of page 5. m is a six-tuple of non-negative intgers and n is a pos. integer. Try
#Tmn([0,1,2,3,4,5],5);
Tmn:=proc(m,n) local i:
binomial(3*n+m[1],m[1])*mul(binomial(2*n,m[i]),i=2..nops(m)):
end:

#Pn(n): the set of primes defined in Lemma2. Try
#Pn(11);
Pn:=proc(n) local p,gu,lu:
gu:={}:
p:=nextprime(max(5,trunc(sqrt(3*n))+1 )):

while p<=2*n do

lu:= n/p-trunc(n/p):

if lu>=1/2 and lu<2/3 then
 gu:=gu union {p}:
fi:

p:=nextprime(p):

od:

gu:
 
end:

#IsGoodpj(p,j,n):Is p divisible by T(m) for all m in M_j?Try:
#IsGoodpjn(5,3,6);
IsGoodpjn:=proc(p,j,n) local gu,gu1:
gu:=Mjn(j,n):

for gu1 in gu do
 if not type(Tmn(gu1,n)/p,integer) then
  print(gu1 , `is a counterexample `):
  RETURN(false):
 fi:
od:

true:

end:

#OrdpTmn(p,m,n): The lower bound for the p-adic order of Tmn(m,n). Try:
#OrdpTmn(11,[0,4,5,3,2,1],30);
OrdpTmn:=proc(p,m,n)  local i:
if not (type(p,prime) and type(m,list) and nops(m)=6 and type(n,integer) and n>0) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

trunc((3*n+m[1])/p)-trunc(3*n/p)-trunc(m[1]/p)+add(trunc(2*n/p)-trunc((2*n-m[i])/p)-trunc(m[i]/p),i=2..6):
end:


#CheckOriL2(n): checks the original Lemma 2. Try:
#CheckOriL2(10);
CheckOriL2:=proc(n) local gu,p,lu,j,lu1:
gu:=Pn(n):

print(`n is`, n):

print(`The set of good primes is`,gu):

for p in gu do

print(`Examining the prime `, p):

print(`Let's take j to be`, p):

j:=p:

lu:=Mjn(j,n) :

print(`Mjn(j,n) has`, nops(lu), `elements , let's see whether the order of p in  each of them is stictly positive `):

for lu1 in lu do

if OrdpTmn(p,lu1,n)<=0 then
 print(lu1, `is a counterexample `):
fi:

od:
od:

end:



#CheckL2(n): checks Lemma 2. Try:
#CheckL2(10);
CheckL2:=proc(n) local gu,p,lu,j:
gu:=Pn(n):


lu:=Alist(n):


for p in gu do

for j from 0 by p to 3*n do
  if not type(lu[j+1]/p, integer) then
   print(`We found a counter-example to Lemma 2, namely (p,j)=`, (p,j) ):
  RETURN(false):
  fi:
od:


od:


lu:=AlistNeg(n):


for p in gu do

for j from  p by p to 4*n do
  if not type(numer(lu[j])/p, integer) then
   print(`We found a counter-example to Lemma 2, namely (p,j)=`, (p,-j) ):
  RETURN(false):
  fi:
od:


od:


true:

end:





#CheckL2v(n): checks Lemma 2. Verbose version. Try:
#CheckL2v(10);
CheckL2v:=proc(n) local gu,p,lu,j:
gu:=Pn(n):

print(`n is`, n):

lu:=Alist(n):

print(`The set of good primes is`,gu):

for p in gu do

print(` p is `, p):

for j from p by p to 3*n do
 print(` j is `,j):
 print(`A[j]/p is `, lu[j+1]/p):
  if not type(lu[j+1]/p, integer) then
   print(`We found a counter-example to Lemma 2, namely (p,j)=`, (p,j) ):
  RETURN(false):
  fi:
od:


od:

true:

end:

#CheckL2G(N): checks Lemma 2 for all n from 1 to N do
#CheckL2G(100);
CheckL2G:=proc(N) local n:

for n from 1 to N do
 if not CheckL2(n) then
  print(`n=`, n, `is a counterexample `):
  RETURN(false):
 fi:
od:
true:
end:




#Alist(n): the list of length 3*n+1 consisting of the A_j from j=0 to 3*n. Try:
#Alist(5);
Alist:=proc(n) local R,u,j:

R:=5*(u-5)^(2*n)*(u^4-20*u^3+156*u^2-560*u+800)^(2*n)/(10-u)^(3*n+1):
R:=taylor(R,u=0,3*n+2):
[seq(coeff(R,u,3*n-j),j=0..3*n)]:
end:



#TayR(n,u): the Taylor expansion of the rational function 5*(u-5)^(2*n)*(u^4-20*u^3+156*u^2-560*u+800)^(2*n)/(10-u)^(3*n+1) up to 3n+1 terms. Try:
#TayR(5,u);
TayR:=proc(n,u) local R:
option remember:
R:=5*(u-5)^(2*n)*(u^4-20*u^3+156*u^2-560*u+800)^(2*n)/(10-u)^(3*n+1):
R:=taylor(R,u=0,3*n+2):
add(coeff(R,u,i)*u^i,i=0..3*n+1):
end:

#Anj(n,j): A_j The coefficient of 1/(x+5)^(j+1) in the partial fration expansion of R(x)
Anj:=proc(n,j)  local u:
coeff(TayR(n,u),u,3*n-j):
end:






#######to be deleted
#TayRmodp(n,u,p): the Taylor expansion of the rational function 5*(u-5)^(2*n)*(u^4-20*u^3+156*u^2-560*u+800)^(2*n)/(10-u)^(3*n+1) up to 3n+1 terms
#modlulo p . In two parts. Try:
#TayRmodp(10,u,13);
#TayRmodp:=proc(n,u,p) local R1,R2,r1,r2,k1,k2:
#2option remember:

#r1:=2*n mod p:
#r2:=3*n mod p:

#2k1:=(2*n-r1)/p:
#k2:=(3*n-r2)/p:

#R1:=5*(u-5)^(k1*p)*(u^4-20*u^3+156*u^2-560*u+800)^(k1*p)/(10-u)^(k2*p):
#R2:=(u-5)^(r1)*(u^4-20*u^3+156*u^2-560*u+800)^(r1)/(10-u)^(r2+1):
#2R1:=taylor(R1,u=0,3*n+2) mod p:
#R2:=taylor(R2,u=0,3*n+2) :

#[R1,R2]:

#end:
#######delete

#Np(p): Inputs a prime p, larger than 5, outputs the set of n such that n<p^2/3 and 1/2<={n/p}<2/3*p. Try:
#Np(17);
Np:=proc(p) local gu, n,lu:

if not isprime(p) and p>5 then
 print(p, `should have been a prime larger than 5`):
 RETURN(FAIL):
fi:

gu:={}:
for n from 1 to trunc(p^2/3) do
  lu:=n/p-trunc(n/p):

  if lu>=1/2 and lu<2/3 then
    gu:=gu union {n}:
  fi:
od:

gu:

end:


#NpG(p,G): Inputs a prime p, larger than 5, outputs the set of n such that n<G and 1/2<={n/p}<2/3*p. Try:
#NpG(17,trunc(17^2/3));
NpG:=proc(p,G) local gu, n,lu:

if not isprime(p) and p>5 then
 print(p, `should have been a prime larger than 5`):
 RETURN(FAIL):
fi:

gu:={}:
for n from 1 to G  do
  lu:=n/p-trunc(n/p):

  if lu>=1/2 and lu<2/3 then
    gu:=gu union {n}:
  fi:
od:

gu:

end:


#CheckL2pGslow(p,G): Inputs a prime p, and a positive integer G, checks the generalized Lemma 2 for all n in NpG(p,G). Try:
#CheckL2pGslow(17,trunc(17^2/3));
CheckL2pGslow:=proc(p,G) local gu,n,j,u:
gu:=NpG(p,G):

for n in gu do
 if {seq(coeff(TayR(n,u) mod p,u,3*n-j*p),j=0..trunc(3*n/p))}<>{0} then
  print(`p=`, p, `and n=`, n, `is a counterexample to the generalized Lemma 2`):
  RETURN(false):
 fi:
od:
true:
end:

#rRangep(p): Inputs a prime p, outouts the set of integers between p/2 and 2*3/p. Try;
#rRangep(11);
rRangep:=proc(p) local i:
{seq(i,i=trunc(p/2)+1..trunc(2*p/3))}:
end:


#TayR1modp(u,p,K): the Taylor expansion of 
#R1=(u^4-20*u^3+156*u^2-560*u+800)^2*(u-5)^2/(10-u)^3 up to u^(K+1) mod p. Try:
#TayR1modp(u,11,30);
TayR1modp:=proc(u,p,K) local R1,i:
option remember:
R1:=u^4-20*u^3+156*u^2-560*u+800 mod p:
R1:=R1^2*(u-5)^2/(10-u)^3 :

R1:=taylor(R1,u=0,K+1) mod p:
add(coeff(R1,u,i)*u^i,i=0..K):
end:


#TayRfastP(u,p,n,K): The first K coefficients of the Taylor exapnsion of
#R1^n where
#R1=(u^4-20*u^3+156*u^2-560*u+800)^2*(u-5)^2/(10-u)^3 up to u^(K+1) mod p. Try:
#TayRfastP(u,7,20,30);
TayRfastP:=proc(u,p,n,K) local mu,gu,i:
option remember:
mu:=TayR1modp(u,p,K):
if n=1 then
 RETURN(mu):
elif n mod 2=0 then
gu:=TayRfastP(u,p,n/2,K)^2:
gu:=taylor(gu,u=0,K+1) mod p:
RETURN(add(coeff(gu,u,i)*u^i,i=0..K)):
else
gu:=TayRfastP(u,p,n-1,K)*mu:
RETURN(add(coeff(gu,u,i)*u^i,i=0..K)):
fi:
end:



#TayRfast(n,u,p): Fast version of the Taylor expansion up to u^(3*n+1)
#5*(u-5)^(2*n)*(u^4-20*u^3+156*u^2-560*u+800)^(2*n)/(10-u)^(3*n+1):
# Try
#TayRfast(20,u,7);
TayRfast:=proc(n,u,p) local gu,mu,i,K:
K:=3*n+1:
gu:=TayRfastP(u,p,n,K):
mu:=taylor(5/(10-u),u,K+1) mod p:
mu:=add(coeff(mu,u,i)*u^i,i=0..K):
gu:=expand(gu*mu) mod p:
gu:=add(coeff(gu,u,i)*u^i,i=0..K):
end:




#CheckL2pG(p,G): Inputs a prime p, and a positive integer G, checks the generalized Lemma 2 for all n in NpG(p,G). Try:
#CheckL2pG(17,trunc(17^2/3)); Fast version
CheckL2pG:=proc(p,G) local gu,n,j,u:
gu:=NpG(p,G):

for n in gu do
 if {seq(coeff(TayRfast(n,u,p),u,3*n-j*p),j=0..trunc(3*n/p))}<>{0} then
  print(`p=`, p, `and n=`, n, `is a counterexample to the generalized Lemma 2`):
  RETURN(false):
 fi:
od:
true:
end:



#Anpj(n,p,j): the A_j mod p in the n case
#Anpj(5,9,0);
Anpj:=proc(n,p,j) local u,R1,R2,i:
option remember:

R1:=expand((u-5)^(2)*(u^4-20*u^3+156*u^2-560*u+800)^(2)/u^3) mod p:
R2:=5/(10-u)^(3*n+1):
R1:=expand(R1^n):

R2:=taylor(R2,u=0,3*n+1):
R2:=add(coeff(R2,u,i)*u^i,i=0..3*n):
coeff(expand(R1*R2),u,-j) mod p:

end:





#CheckL2AOld(n,N): checks the statement in p. 7, lines 4-7 of
# http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/pimeas.pdf for
#n, all primes in Pn(n) and 0<=n0,n1,n2<=N. Try:
#CheckL2(10,10);
CheckL2AOld:=proc(n,N) local gu,n0,n1,n2,p,j:
gu:=Pn(n):

for p in gu do
 for j from 0 by p to 3*n do
   for n0 from 0 to N do
   for n1 from 0 to N do
   for n2 from 0 to N do
    
 if binomial(2*n,n1)*binomial(2*n,n2)*binomial(2*n+2*n1+2*n2,n0)*binomial(3*n+3*n-j-n0,3*n) mod p<>0 then
  print(`(n,p,j,n0,n1,n2)=`, n,p,j,n0,n1,n2 , `is a counterexample `):
print(`Indeed  binomial(2*n,n1)*binomial(2*n,n2)*binomial(2*n+2*n1+2*n2,n0)*binomial(3*n+3*n-j-n0,3*n) equals then`):
print( binomial(2*n,n1)*binomial(2*n,n2)*binomial(2*n+2*n1+2*n2,n0)*binomial(3*n+3*n-j-n0,3*n)):
print(`and this is not divisible by`, p):
  RETURN(false):
 fi:
  od:
  od:
  od:
 od:
od:
true:

end:



B:=proc(n) local u: 
option remember:
coeff(taylor(5*(u-5)^(2*n)*(u^4-20*u^3+156*u^2-560*u+800)^(2*n)/(10-u)^(3*n+1),u=0,3*n+1),u,3*n):end:

Bf:=proc(n1) local N,ope,n,i:
option remember:

ope:=
-2*(2*n-1)*(440102548665*n^8-3486847561002*n^7+11500210953774*n^6-20446104406268*n^5+21171628331761*n^4-12833781780666*n^3+4314166905128*n^2-697150445344*n+
37917022080)/(4*n-1)/(3*n-1)/(3*n-2)/(4*n-3)/n/(559455*n^4-3313544*n^3+7262717*n^2-6977452*n+2478688)/N-64/3*(2*n-1)*(2*n-3)*(n-1)*(1208982255*n^6-7160568584*n^5+
15750386413*n^4-15959336682*n^3+7602260438*n^2-1520500204*n+91752528)/(4*n-1)/(3*n-1)/(3*n-2)/(4*n-3)/n/(559455*n^4-3313544*n^3+7262717*n^2-6977452*n+2478688)/N^2-\
1024/3*(2*n-1)*(2*n-3)*(2*n-5)*(n-1)*(n-2)*(559455*n^4-1075724*n^3+678815*n^2-154830*n+9864)/(4*n-1)/(3*n-1)/(3*n-2)/(4*n-3)/n/(559455*n^4-3313544*n^3+7262717*n^2-\
6977452*n+2478688)/N^3:

if n1<=2 then
 RETURN(B(n1)):
else
add(subs(n=n1,coeff(ope,N,-i))*Bf(n1-i),i=1..3):
fi:

end:

#Checkj0(N):checks Lemma2 for j=0 and all n<=N. Try:
#Checkj0(1000);
Checkj0:=proc(N) local i,p:

evalb({seq(type(Bf(i)/mul(p,p in Pn(i)),integer),i=1..N)}={true}):

end:





#Rnx(n,x): The rational function R_n(x) defined in the first line Equation (2) of the ZZ paper. Try:
#Rnx(3,x);
Rnx:=proc(n,x):
factor(evalc(
5*x^(2*n)*(x+1+2*I)^(2*n)*(x+1-2*I)^(2*n)*(x-1+2*I)^(2*n)*(x-1-2*I)^(2*n)/(5+x)^(3*n+1)/(5-x)^(3*n+1)
)):
end:


#RnxS(n,x): The rational function R_n(x) defined in the second line Equation (2) of the ZZ paper, not using complex numbers. Try:
#RnxS(3,x);
RnxS:=proc(n,x):
normal(evalc(
5*x^(2*n)*
(x^4+6*x^2+25)^(2*n)
/(5+x)^(3*n+1)/(5-x)^(3*n+1)
)):
end:

#CheckEq2(N): Checks that the two expressions for R_n(x) given in Equation (2) are correct for n from 1 to N do
#CheckEq2(20);
CheckEq2:=proc(N) local n1,x:
evalb({seq(normal(Rnx(n1,x)-RnxS(n1,x)),n1=1..N)}={0}):
end:




#CheckEq5c(n): CheckEq5b: cheks the first line of Eq(5) for n. Try:
#CheckEq5c(4);
CheckEq5c:=proc(n) local gu, j,x:

gu:=Rnx(n,x):
gu:=normal((x+5)^(3*n+1)*gu):
evalb([seq(subs(x=-5,diff(gu,x$(3*n-j))/(3*n-j)!  ),j=0..3*n-1),subs(x=-5,gu)]=Alist(n)):
end:



#CheckEq5b(n): CheckEq5b(n): checks the second line of equation (5) of the ZZ paper for n. Try:
#CheckEq5b(4);
CheckEq5b:=proc(n) local gu,Aj,m,ku,i1,j,lu:
gu:=Alist(n):

for j from 0 to 3*n do
 Aj:=gu[j+1]:
 lu:=Mjn(j,n):
 ku:= 5*add(WtEq5(m,n),   m in lu):
 ku:=simplify(evalc(ku)):
 if ku<>Aj then
  print(`When n=`, n, `and j =`, j, `Equation (5)  (second line) is false `):
  print(`Aj is`, Aj, `but the right side is`, ku):
  print(`The ratio is`, Aj/ku):
  RETURN(false):
 fi:
od:

true:

end:

#CheckEq5a(n): CheckEq5a(n): checks the penultimate line of equation (5) of the ZZ paper for n. Try:
#CheckEq5a(4);
CheckEq5a:=proc(n) local gu,Aj,m,ku,i1,j,lu:
gu:=Alist(n):

for j from 0 to 3*n do
 Aj:=gu[j+1]:
 lu:=Mjn(j,n):
 ku:= 5*add((-1)^(add(m[i1],i1=2..6))*Tmn(m,n)*
10^(-3*n-1-m[1])*5^(2*n-m[2])*(4-2*I)^(2*n-m[3])*(4+2*I)^(2*n-m[4])*(6-2*I)^(2*n-m[5])*(6+2*I)^(2*n-m[6]),
   m in lu):
 ku:=simplify(evalc(ku)):
 if ku<>Aj then
  print(`When n=`, n, `and j =`, j, `Equation (5)  (penultimate line) is false `):
  print(`Aj is`, Aj, `but the right side is`, ku):
  print(`The ratio is`, Aj/ku):
  RETURN(false):
 fi:
od:

true:

end:

#CheckEq5Old(n): CheckEq5(n): checks equation (5) of the ZZ paper for n. Try:
#{seq(CheckEq5Old(n1),n1=1..10)}; Try:
#CheckEq5Old(4);
CheckEq5Old:=proc(n) local gu,Aj,m,ku,i1,j,lu:
gu:=Alist(n):

for j from 0 to 3*n do
 Aj:=gu[j+1]:
 lu:=Mjn(j,n):
 ku:= add((-1)^(add(m[i1],i1=2..6))*Tmn(m,n)*2^(4*n-1+j+m[2])*(1-I)^(-m[5])*(1+I)^(-m[6])*5^j*(2+I)^(m[3]+m[6])*(2-I)^(m[4]+m[5]), m in lu):
 ku:=simplify(evalc(ku)):
 if not member(ku/Aj,{-1,1}) then
  print(`When n=`, n, `and j =`, j, `Equation (5) is false `):
  print(`Aj is`, Aj, `but the right side is`, ku):
    print(`The ratio is`, Aj/ku):
  RETURN(false):
 fi:
od:

true:

end:



#CheckEq5(n): CheckEq5(n): checks equation (5) of the ZZ paper for n. Try:
#{seq(CheckEq5(n1),n1=1..10)}; Try:
#CheckEq5(4);
CheckEq5:=proc(n) local j:

evalb([seq(Anj5(n,j),j=0..3*n)]=Alist(n)) and evalb([seq(Anj5(n,-j),j=1..4*n-1)]=AlistNeg(n)):



end:




#WtEq5S(m,n): Given  symbols m and n,outputs the symbolic expression for the summand of the second line of Eq. (5)
#Try:
#WtEq5S(m,n):
WtEq5S:=proc(m,n) local x:
5*
subs(x=-5,diff((5-x)^(-3*n-1),x$m[0]))/m[0]!*
subs(x=-5,diff(x^(2*n),x$m[1]))/m[1]!*
subs(x=-5,diff((x+1+2*I)^(2*n),x$m[2]))/m[2]!*
subs(x=-5,diff((x+1-2*I)^(2*n),x$m[3]))/m[3]!*
subs(x=-5,diff((x-1+2*I)^(2*n),x$m[4]))/m[4]!*
subs(x=-5,diff((x-1-2*I)^(2*n),x$m[5]))/m[5]!:
end:


diff1:=proc(f,x,i):
if i=0 then
 f:
else
diff(f,x$i):
fi:
end:

#WtEq5(m,n): Given a six-tuple  and n,outputs the symbolic expression for the summand of the second line of Eq. (5)
#Try:
#WtEq5(m,n):
WtEq5:=proc(m,n) local x:
subs(x=-5,diff1((5-x)^(-3*n-1),x,m[1]))/m[1]!*
subs(x=-5,diff1(x^(2*n),x,m[2]))/m[2]!*
subs(x=-5,diff1((x+1+2*I)^(2*n),x,m[3]))/m[3]!*
subs(x=-5,diff1((x+1-2*I)^(2*n),x,m[4]))/m[4]!*
subs(x=-5,diff1((x-1+2*I)^(2*n),x,m[5]))/m[5]!*
subs(x=-5,diff1((x-1-2*I)^(2*n),x,m[6]))/m[6]!:
end:



#CompsG(N,A): Inputs a positive ineger N, a list  A, or length r, say ,of non-negative integers outputs
#all the vectors of length r [v1, ..., vr] such that v[i]<=A[i] and that add-up to N
#CompsG(5,[2,2,3]);
CompsG:=proc(N,A) local r,gu,i,lu,lu1:
option remember:

if not type(A,list) then
 print(A, `must be a list`):
 RETURN(FAIL):
fi:

if N<0 then
 RETURN({}):
fi:

r:=nops(A):

if N=0 then
 RETURN({[0$r]}):
fi:



if r=1  then
  if N<=A[1] then
   RETURN({[N]}):
  else
   RETURN({}):
  fi:
fi:

gu:={}:

for i from 0 to A[r] do
 lu:=CompsG(N-i,[op(1..r-1,A)]):
 gu:=gu union {seq([op(lu1),i],lu1 in lu)}:
od:

gu:

end:

#CheckEq5A(n): checks the last line of p.5. Try:
#Check5A(4);
CheckEq5A:=proc(n) local j,m4,m5,gu:

for j from 0 to 3*n do
 for m4 from 0 to min(3*n-j,2*n) do
  for m5 from 0 to min(3*n-j-m4,2*n) do
  gu:=2^ceil((3*n-j)/2)*(1-I)^(-m4)*(1+I)^(-m5):
  gu:=evalc(gu):
  if not (type(coeff(gu,I,0),integer) and type(coeff(gu,I,1),integer)) then
    print(`(n,j,m4,m5)=`, (n,j,m4,m5), `is a counterexample `):
    RETURN(false):
  fi:

 od:
od:
od:

true:

end:



#CheckEq5B(n): checks the first line of p.6. Try:
#Check5B(4);
CheckEq5B:=proc(n) local j,gu,m,lu:

for j from 0 to 3*n do
gu:=Mjn(j,n):

 for m in gu do
   lu:=2^(-floor((5*n+3*j)/2)+1)*2^(4*n-1+j+m[2])*(1-I)^(-m[5])*(1+I)*(-m[6]):
   lu:=evalc(lu):
 if not (type(coeff(lu,I,0),integer) and type(coeff(lu,I,1),integer)) then
    print(`(n,j,m)=`, (n,j,m), `is a counterexample `):
  RETURN(false):
 fi:
od:
od:

true:

end:

 


#Anj5(n,j): The formula for A_j given in Eq(5)
#Try:
#Anj5(3,0);
Anj5:=proc(n,j) local lu,ku,i1,m:
 lu:=Mjn(j,n):

 ku:= add((-1)^(add(m[i1],i1=2..6))*Tmn(m,n)*2^(4*n-1+j+m[2])*(1-I)^(-m[5])*(1+I)^(-m[6])*5^j*(2+I)^(m[3]+m[6])*(2-I)^(m[4]+m[5]), m in lu):
simplify(evalc(ku)):


end:


#Pnx(n,x): The polynomial part of R_n(x). Try:
#Pnx(10,x);
Pnx:=proc(n,x) local gu,mu,i:

mu:=convert(Rnx(n,x),parfrac):
gu:=0:

for  i from 1 to nops(mu) do
 if degree(denom(op(i,mu)),x)=0 then
  gu:=gu+op(i,mu):
 fi:
od:

expand(gu):
end:



#CheckEq6Old(n): Checks Eq. (6) of the ZZ paper. Try:
#CheckEq6Ood(3);
CheckEq6Old:=proc(n) local P,j,k,x:

P:=add((Anj5(n,-k)-add(binomial(j+k-1,j)*Anj5(n,j)/10^(j+k),j=0..3*n))*(x+5)^(k-1),k=1..4*n-1):
evalb(expand(P-Pnx(n,x))=0):
end:





#AlistNegPre(n): The negative part of the A list A[-1], .... A[-4*n]. Try:
#AlistNegPre(3);
AlistNegPre:=proc(n) local gu,x,i:
gu:=Pnx(n,x):
gu:=expand(subs(x=x-5,gu)):
[seq(coeff(gu,x,i-1),i=1..4*n-1)]:
end:


AlistNeg:=proc(n) local gu,j,mu,k:
gu:=AlistNegPre(n):
mu:=Alist(n):

[seq(gu[k]+add(binomial(j+k-1,j)*mu[j+1]/10^(j+k),j=0..3*n),k=1..4*n-1)]:

end:


#Znmj(n,m,j): The discrete function Z(n,m,j) defined in line 8 of p. 7, Try:
#Znmj(4,2,-1);
Znmj:=proc(n,m,j) local n0:
if j>=0 then
add((-2)^n0*binomial(2*n+2*m,n0)*binomial(6*n-j-n0,3*n),n0=0..3*n-j):
else
add((-2)^n0*binomial(2*n+2*m,n0)*binomial(6*n-j-n0,3*n),n0=0..min(2*n+2*m,3*n-j)):
fi:

end:


#Z1nmj(n,m,j): The expression for Z(n,m,j) at the penultimate line of p. 7. Try:
#Z1nmj(5,3,1);
Z1nmj:=proc(n,m,j) local n0:
add((-1)^n0*binomial(2*n+2*m,n0)*binomial(6*n-2*(n+m)-j,3*n-j-n0),n0=0..3*n-j):
end:


#Z2nmj(n,m,j): The expression for Z(n,m,j) in the bottom line of p. 7. Try:
#Z2nmj(5,3,1);
Z2nmj:=proc(n,m,j) local k:
(-1)^(n+m)*add((-1)^(k)*binomial(2*n+2*m,n+m+k)*binomial(4*n-2*m-j,2*n-m-k),k=-(n+m)..(n+m)):
end:




#CheckEq6A(n): checks the correctness of the expression on line 4 for n and all j in {-4*n,..., 3*n}. Try:
#CheckEq6A(10);
CheckEq6A:=proc(n) local j,lu,ku,n1,n2:


lu:=Alist(n):

for j from 0 to 3*n do
ku:=add(add((3-4*I)^(2*n-n1)*(3+4*I)^(2*n-n2)*5^(2*n+2*n1+2*n2+1)*10^(-(6*n-j)-1)*binomial(2*n,n1)*binomial(2*n,n2)*Znmj(n,n1+n2,j),
n2=0..2*n),n1=0..2*n):
 ku:=evalc(ku):
 if ku<>lu[j+1] then
  print(`(n,j)=`, [n,j], `is a counterexample,the ratio being`):
  print(ku/lu[j+1]):
  RETURN(false):
 fi:
od:





lu:=AlistNeg(n):

for j from -1  by -1 to -(4*n-1) do
ku:=add(add((3-4*I)^(2*n-n1)*(3+4*I)^(2*n-n2)*5^(2*n+2*n1+2*n2+1)*10^(-(6*n-j)-1)*binomial(2*n,n1)*binomial(2*n,n2)*Znmj(n,n1+n2,j),
n2=0..2*n),n1=0..2*n):
 ku:=evalc(ku):
 if ku<>lu[-j] then
 print(`(n,j)=`, [n,j], `is a counterexample,the ratio being`):
 print(ku/lu[-j]):
 RETURN(false):
fi:
od:

true:

end:

 





#CheckL2Aold(n,M): checks the statement in the middle  of page 7 about the divisibility of That(n,n1,n2)=binomial(2*n,n1)*binomial(2*n,n2)*Znmj(n,m1+m2,j)
#for all j from 0 to 3*n and all n1,n2<=M. Try:
#CheckL2Aold(4,10);
#CheckL2Aold(10,3);
CheckL2Aold:=proc(n,M) local gu,p,j,n1,n2:

gu:=Pn(n):

for p in gu do
 for j from 0 by p to 3*n do
  for n1 from 0 to M do
   for n2 from 0 to M do
    if binomial(2*n,n1)*binomial(2*n,n2)*Znmj(n,n1+n2,j) mod p<>0 then

     print(`(n,n1,n2,p,j)=`, (n,n1,n2,p,j), `is a counterexample `):
     RETURN(false):
    fi:
   od:
  od:
 od:
od:
    
true:

end:




#CheckL2A(n): checks the statement in the middle  of page 7 about the divisibility of That(n,n1,n2)=binomial(2*n,n1)*binomial(2*n,n2)*Znmj(n,m1+m2,j)
#for all j from 0 to 3*n and all n1,n2<=2*n. Try:
#CheckL2A(10);
#CheckL2A(10);
CheckL2A:=proc(n) local gu,p,j,n1,n2:

gu:=Pn(n):

for p in gu do
 for j from 0 by p to 3*n do
  for n1 from 0 to 2*n do
   for n2 from 0 to 2*n do
    if binomial(2*n,n1)*binomial(2*n,n2)*Znmj(n,n1+n2,j) mod p<>0 then

     print(`(n,n1,n2,p,j)=`, (n,n1,n2,p,j), `is a counterexample `):
     RETURN(false):
    fi:

   od:
  od:
 od:
od:
    
true:

end:





#ordp(n,p): the highest exponent i such that  n/p^i is an integer. Try:
#ordp(250,5);
ordp:=proc(n,p) local i:
for i from 0 while n mod p^i=0  do od:
i:=i-1:
end:


#ordpG(n,p): For a rational number n, the highest exponent i such that  n/p^i is an integer. Try:
#ordpG(19/5,5);
ordpG:=proc(n,p) local i:
ordp(numer(n),p)-ordp(denom(n),p):
end:


#CheckEq6B(N,K): checks that the p-order of N! is floor(N/p) if p>sqrt(N) up to K primes. Try:
#heckEq6B(10,10);
CheckEq6B:=proc(N,K) local p,gu,i:

p:=nextprime(floor(evalf(sqrt(N)))):

for i from 1 to K do
 gu:=ordp(N!,p):
 if gu<>floor(N/p) then
  print(`N=`, N, `p= `, p, `is a counterexample `):
  RETURN(false):
 fi:
od:
true:
end:




#CheckSTAM(N,K): Checks Equation (7) in the ZZ paper for all n<=N, n1<=2*n,  and the first K primes. Try:
#CheckSTAM(10,10);
CheckSTAM:=proc(N,K) local n,n1,p,k:

for n from 1 to N do
for n1 from 1 to 2*n do
 for k from 1 to K do
    p:=ithprime(k):
   if p>evalf(sqrt(2*n)) then

  if ordp(binomial(2*n,n1),p)<>frac(n1/p)+frac((2*n-n1)/p)-frac(2*n/p) then
   print(` (n,n1,p)=`, (n,n1,p), `is a counterexample `):
   RETURN(false):
  fi:
  fi:
 od:
od:
od:
true:

end:

#CheckEq7(N,K): Checks Equation (7) in the ZZ paper for all n<=N, n1<=2*n,  and the first K primes. Try:
#CheckEq7(10,10);
CheckEq7:=proc(N,K) local n,n1,p,w,w1,k:

for n from 1 to N do
for n1 from 1 to 2*n do
 for k from 1 to K do
    p:=ithprime(k):
   if p>evalf(sqrt(2*n)) then
    w:=frac(n/p):

    w1:=frac(n1/p):


  if ordp(binomial(2*n,n1),p)<>floor(2*w)-floor(2*w-w1) then
   print(`When n=`, n):
   print(`When n1=`, n1):
   print(`When p=`, p):
    print(`binomial(2*n,n1) is`,  binomial(2*n,n1)):
    print(`its p-order is`,  ordp(binomial(2*n,n1),p) ):
    print(`The fractional part of n/p, w,  is`, w):
    print(`The fractional part of n1/p, w1,  is`, w1):
     print(`2*w is`, 2*w):
    print(`2*w-w1 is`, 2*w-w1):
    print(`floor(2*w)-floor(2*w-w1) is`, floor(2*w)-floor(2*w-w1)):
   print(`Hence (n,n1,p)=`, (n,n1,p), `is a counterexample `):
   RETURN(false):
  fi:
  fi:
 od:
od:
od:
true:

end:







#PHIn(n):PHI_n defined in the statement of Lemma 3 on p. 9. Try:
#PHIn(10);
PHIn:=proc(n) local gu,p:
gu:=Pn(n):
mul(p, p in gu):
end:

#Ln(n):L_n defined in the statement of Lemma 3 on p. 9. Try:
#Ln(10);
Ln:=proc(n) local i:
lcm(seq(i,i=1..4*n))/PHIn(n):
end:


CheckL3:=proc(N) local n, guP,j:

for n from 1 to N do
guP:=Alist(n):

if not type(guP[1]/PHIn(n), integer) then
 print(`(n,j)=`, (n,0), ` is a counterexample `):
 RETURN(false):
fi:

for j from 1 to 3*n do
 if not type(Ln(n)*guP[j+1]*10*(-j)/j, integer) then
  print(`(n,j)=`, (n,j), ` is a counterexample `):
  RETURN(false):
 fi:
od:



guP:=AlistNeg(n):

for j from 1 to 4*n-1 do
 if not type(Ln(n)*guP[j]*10*(j)/(-j), integer) then
  print(`(n,j)=`, (n,-j), ` is a counterexample `):
  RETURN(false):
 fi:
od:

od:

true:
end:



#CheckL4(N): cheks Lemma 4 on p. 10 for n between 1 and N. Try:
#CheckL4(10);
CheckL4:=proc(N) local n,gu,x,k,Bk,lu:

for n from 1 to N do
gu:=Pnx(n,x):
gu:=expand(evalc(subs(x=x-1-2*I,gu))):

for k from 0 to 4*n-2 do
   Bk:=coeff(gu,x,k):
   lu:=2^(-floor(5*n/2)+ceil(3*k/2)+3)*Bk:
   if not (type(coeff(lu,I,0),integer) and type(coeff(lu,I,1),integer)) then
    print(`(n,k)`, n,k, `2^(-floor(5*n/2)+ floor(3*k/2)+3) * B_k is `, lu):
    print(`Hence (n,k)`, n,k, `is a counterexample to Lemma 4`):
    RETURN(false):
  fi:
od:

od:

true:

end:


#CheckL5(N): cheks Lemma 5 on p. 10 for n from 1 to N. Try:
#CheckL5(10);
CheckL5:=proc(N) local n,gu,x,lu:

for n from 1 to N do
gu:=Pnx(n,x):

gu:=2^(-floor(5*n/2))*Ln(n)*I*int(gu,x=-1-2*I..-1+2*I):
gu:=simplify(evalc(gu)):

if not type(gu, integer) then
 print(`n=`, n, `is a counterexample to Lemma 5`):
 RETURN(false)
fi:

od:

true:

end:


#CheckL6(N): checks Lemma 6 on p. 11 for n from 1 to N. Try:
#CheckL6(10);
CheckL6:=proc(N) local n,gu,x,j,lu:

for n from 1 to N do

lu:=Alist(n):

gu:=add(lu[j+1]*(1/(5+x)^(j+1)+1/(5-x)^(j+1)),j=1..3*n):
gu:=2^(-trunc(5*n/2)+1)*Ln(n)*I*int(gu,x=-1-2*I..-1+2*I):
gu:=simplify(evalc(gu)):
if not type(gu, integer) then
 print(`The value is`, gu ):
 print(`n=`, n, `is a counterexample to Lemma 5`):
 RETURN(false)
fi:

od:

true:

end:

 

#CheckP1(N): Checks Proposition 1 on p.12 from 1 to N. Try:
#CheckP1(20);
CheckP1:=proc(N) local n,gu,x:

for n from 1 to N do
gu:=Rnx(n,x):
gu:=I*(-1)^(n+1)*int(gu,x=-1-2*I..-1+2*I):
gu:=simplify(evalc(gu)):
gu:=2^(-floor(5*n/2)+2)*Ln(n)*gu:

if not (type(coeff(gu,Pi,0),integer) and type(coeff(gu,Pi,1),integer)) then
 print(`when n=`, n, `the value is`, gu):
 print(`Hence this is a counterexample `):
 RETURN(false):
fi:

od:

true:

end:


WAD:=proc(a,b)
if a=0 and b=0 then
 1:
elif a<>0 and b<>0 then
a/b:
else
FAIL:
fi:
end:


#CheckEq7A(N): Checks the alternative expressions for Z(n,m,j) given in the last two lines of p. 7 for
#n,m<=N and j between 0 and 3*n. Try:
#CheckEq7A(10);
CheckEq7A:=proc(N) local j,m,n:

evalb({seq(seq(seq(Znmj(n,m,j)-Z1nmj(n,m,j),j=0..3*n),n=0..N),m=0..N)}={0} and {seq(seq(seq(Znmj(n,m,j)-Z2nmj(n,m,j),j=0..4*n-2*m),n=0..N),m=0..N)}={0}):
end:



#CheckEq7B(K): Checks  the identity in the third line of  p. 8 for 0<=M<=N<=K. Try
#CheckEq7B(10);
CheckEq7B:=proc(K) local M,N,k:evalb({seq(seq(add((-1)^k*binomial(2*N,N+k)*binomial(2*M,M+k),k=-N..N)-(2*N)!*(2*M)!/(N!*(N+M)!*M!),N=0..M),M=0..K)}={0}):end:



#CheckEq7C(K): Checks  the sum in the 5th line and is indeed the coefficient of t^(2*N) in the expression in the 5th line. Try:
#CheckEq7C(10);
CheckEq7C:=proc(K) local M,N,k,j,t:
evalb({
seq(  seq(    seq(         add((-1)^k*binomial(2*N,N+k)*binomial(2*M-j,M+k),k=-N..N)-
(-1)^N*(2*N)!*(2*M-j)!/((N+M)!*(N+M-j)!)*coeff((1+t)^(N+M)*(1-t)^(N+M-j),t,2*N),j=0..2*M),M=0..N),N=0..K)}={0}):
end:





#CheckEq7Dold(N,K): Checks the statement about the p-adic order of (2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)!
#for all0<=n1<=n2<=2*n and  j between 0 and min(4*n-2*n1-2*n2,3*n). Try:
#CheckEq7Dold(10,10);
CheckEq7Dold:=proc(N,K) local n,n1,n2,j, p,w,w1,w2, k,gu,mu,gu1:

for n from 1 to N do
for n1 from 0 to 2*n do
for n2 from 0 to 2*n do
for j from 0 to min(4*n-2*n1-2*n2,3*n) do
 for k from 1 to K do
    p:=ithprime(k):

 gu:=(2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)!:

   gu1:=ordpG(gu,p):
  

    w:=frac(n/p):
    w1:=frac(n1/p):
     w2:=frac(n2/p):
    
    mu:=floor(2*w+2*w1+2*w2)+ floor(4*w-2*w1-2*w2-j/p)-floor(3*w)-floor(3*w-j/p):
   if gu1<>mu then
    print(`(n,n1,n2,p,j)=`, (n,n1,n2,p,j) , `is a counterexample `):
     print(`(2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)! is`, gu):
  print(`The order is`, gu1, `but the expression is`, mu):
    RETURN(false):
   fi:
 od:
od:
od:
od:
od:

true:

end:





#CheckEq7D(N): Checks the statement about the p-adic order of (2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)!
#for all0<=n1<=n2<=2*n and  j between 0 and min(4*n-2*n1-2*n2,3*n). Try:
#CheckEq7D(10,10);
CheckEq7D:=proc(N) local vu,n,n1,n2,j, p,w,w1,w2, gu,mu,gu1:

for n from 1 to N do
vu:=Pn(n):
for n1 from 0 to 2*n do
for n2 from 0 to 2*n do
     for  p in vu do
for j from  0 to min(3*n, 4*n-2*n1-2*n2) do

 gu:=(2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)!:

   gu1:=ordpG(gu,p):
  

    w:=frac(n/p):
    w1:=frac(n1/p):
     w2:=frac(n2/p):

if 2*w-w1>=1 and 2*w-w2>=1 then    
    mu:=floor(2*w+2*w1+2*w2)+ floor(4*w-2*w1-2*w2-j/p)-floor(3*w)-floor(3*w-j/p):
   if gu1<>mu then
    print(`(n,n1,n2,p,j)=`, (n,n1,n2,p,j) , `is a counterexample `):
     print(`(2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)! is`, gu):
  print(`The order is`, gu1, `but the expression is`, mu):
    RETURN(false):
   fi:
 fi:
 od:
od:
od:
od:
od:

true:

end:




#CheckEq7D1(N): Checks the statement about the p-adic order of (2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)!
#for all0<=n1<=n2<=2*n and  j between 0 and min(4*n-2*n1-2*n2,3*n). Try:
#CheckEq7D1(10,10);
CheckEq7D1:=proc(N) local vu,n,n1,n2,j, p,w,w1,w2, gu,mu,gu1:

for n from 1 to N do
vu:=Pn(n):
for n1 from 0 to 2*n do
for n2 from 0 to 2*n do
     for  p in vu do
for j from  0 by p to min(3*n, 4*n-2*n1-2*n2) do

 gu:=(2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)!:

   gu1:=ordpG(gu,p):
  

    w:=frac(n/p):
    w1:=frac(n1/p):
     w2:=frac(n2/p):

if 2*w-w1>=1 and 2*w-w2>=1 then    
    mu:=floor(2*w+2*w1+2*w2)+ floor(4*w-2*w1-2*w2)-2*floor(3*w):
   if gu1<>mu then
    print(`(n,n1,n2,p,j)=`, (n,n1,n2,p,j) , `is a counterexample `):
     lprint(`(2*n+2*n1+2*n2)!*(4*n-2*n1-2*n2-j)!/(3*n)!/(3*n-j)! is`, gu):
     print(`(w,w1,w2)=`, (w,w1,w2)):
    
  print(`The order is`, gu1):
   print(`floor(2*w+2*w1+2*w2)+ floor(4*w-2*w1-2*w2)-2*floor(3*w) is`, mu):
 print(`These are not the same`):

    RETURN(false):
   fi:
 fi:
 od:
od:
od:
od:
od:

true:

end:





#CheckL2Aslow(n): checks the statement in the middle  of page 7 about the divisibility of That(n,n1,n2)=binomial(2*n,n1)*binomial(2*n,n2)*Znmj(n,m1+m2,j)
#for all j from 0 to 3*n and all n1,n2<=2*n. Try:
#CheckL2Aslow(10);
#CheckL2Aslow(10);
CheckL2Aslow:=proc(n) local gu,p,j,n1,n2:

gu:=Pn(n):

for p in gu do
 for j from -(4*n-1)   to 3*n do
 if j mod p=0 then
  for n1 from 0 to 2*n do
   for n2 from 0 to 2*n do
    if binomial(2*n,n1)*binomial(2*n,n2)*Znmj(n,n1+n2,j) mod p<>0 then

     print(`(n,n1,n2,p,j)=`, (n,n1,n2,p,j), `is a counterexample `):
     RETURN(false):
    fi:
    od:
  od:
fi:

 od:
od:
    
true:

end:

#Refusinks11(n): inputs a positive integer n, and  outputs the set of pairs [p,n1], where 0<=n1<=2*n
# such that binomial(2*n,n1) is NOT divisible by p.
#Try Refusinks11(10);
Refusinks11:=proc(n) local gu,mu,n1,p,gu1:
mu:=Pn(n):

gu:=[]:

for p in mu do
 gu1:=[]:
 for n1 from 0 to n do
  if binomial(2*n,n1) mod p<>0 then
    gu1:=[op(gu1),n1]:
  fi:
 od:
  gu:=[op(gu),[p,gu1]]:
od:

gu:

end:




#Refusinks1(n): inputs a positive integer n, and  outputs the set of pairs [p,Listp], such that 
#binomial(2*n,n1) mod p is not divisibile by p
# Refusinks1(10);

Refusinks1:=proc(n) local gu,mu1,i1,i,mu:
gu:=Refusinks11(n):

mu:=[]:


for i from 1 to nops(gu) do
mu1:=gu[i][2]:
mu1:=[op(mu1),op(sort([seq(2*n-mu1[i1],i1=1..nops(mu1))]))]:
mu:=[op(mu),[gu[i][1],mu1] ]:
od:

mu:
end:





#RefusinksB(n): inputs a positive integer n, and  outputs the set of pairs [p,m], such that 
#0<=n1<=2*n,  0<=n2<=2*n, and m=n1+n2, and binomial(2*n,n1)*binomial(2*n,n2)
#is not divisible by p. Try:
# Refusinks1(10);

RefusinksB:=proc(n) local gu,mu,i1,i2,gu1,i:
gu:=Refusinks1(n):

mu:=[]:

for i from 1 to nops(gu) do
 gu1:=gu[i][2]:
 gu1:={seq(seq(gu1[i1]+gu1[i2],i1=1..nops(gu1)),i2=1..nops(gu1))}:
  gu1:=sort(convert(gu1,list)):
  mu:=[op(mu),[gu[i][1],gu1]]:
od:

mu:

end:



#RefusinksZ1(n,i): Inputs a positive integer n and output the list of pairs [p,Lisp] where Listp is the list
#of m such that Z(n,m,i*p) is not divisibile by p. Try
#RefusinksZ1(20,0);
RefusinksZ1:=proc(n,i) local mu,p, gu,m,lu:

mu:=Pn(n):

gu:=[]:
for p in mu do
 lu:=[]:
 for m from 0 to 4*n do
   if Znmj(n,m,p*i) mod p<>0 then
    lu:=[op(lu),m]:
   fi:
 od:
gu:=[op(gu),[p,lu]]:
od:

gu:
end:


#RefusinksZ(n): Inputs a positive integer n and output the list of pairs [p,Lisp] where Listp is the list
#of m such that Z(n,m,i*p) is not divisibile by p for at least one i. Try
#RefusinksZ(20);
RefusinksZ:=proc(n) local mu,p, gu,m,lu,i:

mu:=Pn(n):

gu:=[]:
for p in mu do
 lu:=[]:
 for m from 0 to 4*n do
   if {seq(Znmj(n,m,p*i) mod p,i=-trunc((4*n-1)/p)..trunc(3*n/p) ) }<>{0} then
    lu:=[op(lu),m]:
   fi:
 od:
gu:=[op(gu),[p,lu]]:
od:

gu:
end:


#CheckL2A(n): checks the statement in the middle  of page 7 about the divisibility of That(n,n1,n2)=binomial(2*n,n1)*binomial(2*n,n2)*Znmj(n,m1+m2,j)
#for all j from 0 to 3*n and all n1,n2<=2*n. Fast version. Try
#CheckL2A(10);
CheckL2A:=proc(n) local gu1,gu2,lu,i:

gu1:=RefusinksB(n):
gu2:=RefusinksZ(n):

for i from 1 to nops(gu1) do
 lu:=convert(gu1[i][2],set) intersect convert(gu2[i][2],set):
  if lu<>{} then
   print(`For the member of P_n, the prime`, gu1[i][1], `the following m resuse both the binomial part and and the Z(n,m,j) part `):
   print(lu):
   RETURN(false):
 fi:
od:

true:

end:


#CheckL2Apn(p,n): Checks Lemma2A for a specific p and n Try:
#CheckL2Apn(7,9);
CheckL2Apn:=proc(p,n) local gu,i,gu1,gu2,lu,m:

gu:={}:

for i from 0 to 2*n do
 if not binomial(2*n,i) mod p=0 then
  gu:=gu union {i}:
 fi:
od:

gu:={seq(seq(gu1+gu2,gu1 in gu), gu2 in gu)}:

lu:={}:

for m from 0 to 4*n do
 
if {seq(Znmj(n,m,p*i) mod p,i=-trunc((4*n-1)/p)..trunc(3*n/p) ) }<>{0} then
 lu:=lu union {m}:
fi:
od:


lu intersect gu:

end: