######################################################################
##SuperCongruence.txt: Save this file as  SuperCongruence.txt        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  SuperCongruence.txt<Enter>                   #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: 
 
print(`Created:  Feb. 2019`):
print(` This is SuperCongruence.txt `):
print(`It is the Maple package that accompanies the article `):
print(` Polynomial Reduction and Super Congruences `):
print(`by  Qing-Hu Hou, Uan-Ping Mu, and Doron Zeilberger`):
print(`and also available from Zeilberger's website`):
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: Check33a, CheckHamme1, CheckSun1, En, IsDegen, rf `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Check11a, Check11b, Check22, Check33, Check34, CheckHamme, CheckHammeLong, CheckSun, Findam `):
 print(` `):


elif nops([args])=1 and op(1,[args])=Check11a then
print(`Check11a(p,m,r): Inputs an odd prime p, an odd integer m, and a positive integer r,`):
print(`Checks the first part ( Eq. (1.1)) of Conjecture 1.1 Try:`):
print(`Check11a(5,9,2);`):

elif nops([args])=1 and op(1,[args])=Check11b then
print(`Check11b(p,m,r): Inputs an odd prime p, an odd integer m, and a positive integer r,`):
print(`Checks the second part ( Eq. (1.2)) of Conjecture 1.1 Try:`):
print(`Check11b(5,9,2);`):

elif nops([args])=1 and op(1,[args])=Check22 then
print(`Check22(P,x,k): checks Lemma 2.2 of the Hou-Mu-Zeilberger paper for the pair of polynomials P=[a,b] in the variable k.`):
print(`and a polynomial x in k.`):
print(`Try:`):
print(`Check22([k^3,(k+2)^3],k^3+1,k);`):
print(`Check22([k^3,-(k+2)^3],k^5+1,k);`):

elif nops([args])=1 and op(1,[args])=Check33 then
print(`Check33(k,alpha,r,m,d1): checks Theorem 3.3 of the paper with the rational number alpha, and r and m  positive integers.`):
print(`looks for a polynomial x(k) of degree up to Maxd1.`):
print(`Try: `):
print(`Check33(k,2/5,3,2,20);`):

elif nops([args])=1 and op(1,[args])=Check33a then
print(`Check33a(k,alpha,r,m,d1): checks Theorem 3.3 of the paper with the rational number alpha, and r and m  positive integers.`):
print(`looks for a polynomial x(k) of degree d1.`):
print(`Try: `):
print(`Check33a(k,2/5,3,2,2);`):


elif nops([args])=1 and op(1,[args])=Check34 then
print(`Check34(m,NuP,Maxd): inputs an odd positive integer NuP (for the number of primes to check)`):
print(` to  test, and Maxd for the call of Check33 (make it big enough)`):
print(`and outputs pair [a_m,c_m], where a_m, c_m promised by Theorem 3.4 in the paper`):
print(`and checks that indeed the congruence, mod p^4 is right for all primes<=5 and less than MaxP. Try: `):
print(`Check34(7,10,20);`):

elif nops([args])=1 and op(1,[args])=CheckHamme then
print(`CheckHamme(N): checks's Hamme's super-congruence (top of p.2 of the paper) for the first N odd primes. Try:`):
print(`CheckHamme(20);`):

elif nops([args])=1 and op(1,[args])=CheckHammeLong then
print(`CheckHammeLong(N): checks's Hamme's second super-congruence (line 10 of p.2 of the paper) for the first N  primes 1 (mod 4) (starting at p=5). Try`):
print(`CheckHammeLong(30);`):

elif nops([args])=1 and op(1,[args])=CheckHamme1 then
print(`CheckHamme1(p): checks's Hamme's super-congruence (top of p.2 of the paper) for the prime p. Try:`):
print(`CheckHamme1(11);`):

elif nops([args])=1 and op(1,[args])=CheckSun then
print(`CheckSun(N): checks's Sun's super-congruence (line 7 of p.2 of the paper) for the first N odd primes, starting with p=5. Try:`):
print(`CheckSun(20);`):

elif nops([args])=1 and op(1,[args])=En then
print(`En(n): The n-th Eulerian number. Try:`):
print(`En(5);`):

elif nops([args])=1 and op(1,[args])=Findam then
print(` Findam(p,m):Finds the value of a_m in conjecture 1.1 of the Hou-Mu-Zeilberger paper. Note that is also depends on p. Try:`):
print(`Findam(5,3);`):


elif nops([args])=1 and op(1,[args])=IsDegen then
print(`IsDegen(P,k): inputsa pair of  polynomials P=[a,b] of the variable k, and outputs true iff the pair P=[a,b] is degenerated.`):
print(`Try: `):
print(`IsDegen([k^3,(k+1)^2],k);`):

elif nops([args])=1 and op(1,[args])=rf then
print(`rf(a,k): the raising factorial a(a+1) ..(a+k-1). Try:`):
print(`rf(1/2,10);`):

print(``):


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


#rf(a,k): the raising factorial a(a+1) ..(a+k-1). Try:
#rf(1/2,10);
rf:=proc(a,k) local i:
mul(a+i,i=0..k-1):
end:

#CheckHamme1(p): checks's Hamme's super-congruence (top of p.2 of the paper) for the prime p. Try:
#CheckHamme1(11);
CheckHamme1:=proc(p) local k:

if not (type(p,integer) and p>=3 and isprime(p)) then
 print(p, `is not an odd prime `):
 RETURN(FAIL):
fi:

[add((-1)^k*(4*k+1)*(rf(1/2,k)/k!)^3,k=0..(p-1)/2) mod p^3, (-1)^((p-1)/2)*p mod p^3]:

end:

#CheckHamme(N): checks's Hamme's super-congruence (top of p.2 of the paper) for the first N odd primes. Try
#CheckHamme(30);
CheckHamme:=proc(N) local i:

{seq(CheckHamme1(ithprime(i)),i=2..N+1)}:
end:


#En(n): The n-th Eulerian number. Try:
#En(5);
En:=proc(n) local x:
n!*coeff(taylor(2/(exp(x)+exp(-x)),x=0,n+1),x,n):
end:


#CheckSun1(p): checks's Sun's super-congruence (line 7 of p.2 of the paper) for the prime p. Try:
#ChecSun1(11);
CheckSun1:=proc(p) local k:

if not (type(p,integer) and p>=3 and isprime(p)) then
 print(p, `is not an odd prime `):
 RETURN(FAIL):
fi:

evalb(numer(
add((-1)^k*(4*k+1)*(rf(1/2,k)/k!)^3,k=0..(p-1)/2)- (-1)^((p-1)/2)*p -p^3*En(p-3)
) 
mod p^4 =0):

end:

#CheckSun(N): checks's Hamme's super-congruence (top of p.2 of the paper) for the first N odd primes (starting at p=5). Try
#CheckSun(30);
CheckSun:=proc(N) local i:

evalb({seq(CheckSun1(ithprime(i)),i=3..N+2)}={true}):
end:

#CheckHammeLong1(p): checks's Hamme's second super-congruence (line 10 p.2 of the paper) for the prime p. Try:
#CheckHammeLong1(11);
CheckHammeLong1:=proc(p) local k:

if not (type(p,integer) and p>=3 and isprime(p) and p mod 4=1) then
 print(p, `is not an odd prime 1 mod 4`):
 RETURN(FAIL):
fi:

evalb(numer(add((4*k+1)*(rf(1/2,k)/k!)^4,k=0..(p-1)/2)- p) mod p^4=0):

end:

#CheckHammeLong(N): checks's Hamme's second super-congruence (line 10 of p.2 of the paper) for the first N  primes 1 (mod 4) (starting at p=5). Try
#CheckHammeLong(30);
CheckHammeLong:=proc(N) local co,p:

p:=5:
co:=0:

while co<=N do

if p mod 4=1 then

co:=co+1:
 if not CheckHammeLong1(p) then 
   RETURN(false):
 fi:
fi:

p:=nextprime(p):
od:

true:

end:

# Findam(p,m):Checks that  the value of a_m in conjecture 1.1 of the Hou-Mu-Zeilberger paper. Note that is also depends on p. Try
#Findam(5,3);

Findam:=proc(p,m) local am,gu1,gu2,k:
if not ( type((m-1)/2,integer) and isprime(p) and p>2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

am:=Check34(m,3, 30)[1]:

gu1:=add((-1)^k*(4*k+1)^m*(rf(1/2,k)/k!)^3,k=0..(p-1)/2) mod p^3:

gu2:=(-1)^((p-1)/2)*p  mod p^3:

gu2*am-gu1 mod p^3:


end:

#Check11a(p,m,r): Inputs an odd prime p, an odd integer m, and a positive integer r,
#Checks the first part  (Eq. (1.1)) of Conjecture 1.1 Try:
#Check11a(5,9,2);
Check11a:=proc(p,m,r) local am,k:
if not ( type((m-1)/2,integer) and isprime(p) and p>2 and r>1) then
 print(`Bad input`):
 RETURN(FAIL):
fi:


am:=Findam(p,m):

evalb((add((-1)^k*(4*k+1)^m*(rf(1/2,k)/k!)^3,k=0..(p^r-1)/2 ) -am*p^r*(-1)^((p-1)/2*r)  ) mod p^(r+2)=0):

end:


#Check11b(p,m,r): Inputs an odd prime p, an odd integer m, and a positive integer r,
#Checks the second (Eq. (1.2)) of Conjecture 1.1 Try:
#Check11b(5,9,2);
Check11b:=proc(p,m,r) local bm,k:
if not ( type((m-1)/2,integer) and isprime(p) and p>2 and r>1) then
 print(`Bad input`):
 RETURN(FAIL):
fi:


bm:=Findbm(p,m):

evalb((add((4*k+1)^m*(rf(1/2,k)/k!)^4,k=0..(p^r-1)/2 ) -bm*p^r ) mod p^(r+3) = 0):

end:

#IsDegen(P,k): inputsa pair of  polynomials P=[a,b] of the variable k, and outputs true iff the pair P=[a,b] is degenerate.
#Try
#IsDegen([k^3,k^2],k);
IsDegen:=proc(P,k) local a,b,u,m0:
a:=P[1]: b:=P[2]:
u:=expand(a- subs(k=k-1,b)):
m0:=-lcoeff(u,k)/lcoeff(a,k):

if type(m0, integer) and m0>=0 and degree(u,k)=degree(a,k)-1 then
 true:
else
 false:
fi:
end:


#Check22(P,x,k): checks Lemma 2.2 of the Hou-Mu-Zeilberger paper for the pair of polynomials P=[a,b] in the variable k.
#and a polynomial x in k.
#Try:
#Check22([k^3,(k+2)^3],k^3+1,k);
#Check22([k^3,-(k+2)^3],k^5+1,k);
Check22:=proc(Pa,x,k) local a,b,p,u,m0,d:
a:=Pa[1]: b:=Pa[2]:
 u:=expand(a- subs(k=k-1,b)):

 d:=max(degree(u,k), degree(a,k)-1):

p:=expand(a*subs(k=k+1,x)-subs(k=k-1,b)*x):

print(` p is`, p):

if IsDegen(Pa,k) then
 print(`It is a degenerated pair`):
else
 print(`It is not  a degenerated pair`):
fi:

if IsDegen(Pa,k) then
 m0:=-lcoeff(u,k)/lcoeff(a,k):
 if degree(x,k)=m0 then
  print(`p is `,  p):
  print( `d+m0 is`, d+m0):
   RETURN(evalb(degree(p,k)<d+m0)):
 fi:
fi:


    RETURN(evalb(degree(p,k)=d+degree(x,k))) :

end:

#Check33a(k,alpha,r,m,d1): checks Theorem 3.3 of the paper with the rational number alpha, and r and m  positive integers.
#looks for a polynomial x(k) of degree d1.
#Try:
#Check33a(2/5,3,2,2);
Check33a:=proc(k,alpha,r,m,d1) local D1, eq,var,x,i,a,gu,tkRat,X,var1:

if not (type(alpha,fraction) and type(r,integer) and r>=1 and type(m,integer) and m>=1) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

D1:=denom(alpha):


tkRat:=-(alpha+k)^r/(k+1)^r:

X:=add(x[i]*k^i,i=0..d1):

gu:=(2*D1*k+D1*alpha)^m
- add(a[i]*(2*D1*k+D1*alpha)^i,i=0..r-1)+ 2^(r-1)*(D1*k)^r*subs(k=2*D1*k,X)
-2^(r-1)*(D1*(k+1))^r*subs(k=2*D1*(k+1),X)*tkRat:
gu:=numer(expand(normal(gu))):

var:={seq(x[i],i=0..d1), seq(a[i],i=0..r-1)}:
eq:={seq(coeff(gu,k,i),i=0..degree(gu,k))}:

var1:=solve(eq,var):

if var1=NULL then
  RETURN(FAIL):
fi:
[subs(var1,X),subs(var1,[seq(a[i],i=0..r-1)])]:
end:



 


#Check33(k,alpha,r,m,Maxd): checks Theorem 3.3 of the paper with the rational number alpha, and r and m  positive integers.
#looks for a polynomial x(k) of degree up to Maxd.
#Try:
#Check33(k,2/5,3,2,10);
Check33:=proc(k,alpha,r,m,Maxd) local d1,gu,gu2,i1:

for d1 from 0 to Maxd do
  gu:=Check33a(k,alpha,r,m,d1):
  if gu<>FAIL then
   gu2:=gu[2]:
   if m mod 2=0 then
        if nops(gu2)>=2 then
    if {seq(gu2[2*i1],i1=1..trunc(nops(gu2)/2))}<>{0} then
  print(   {seq(gu2[2*i1-1],i1=1..trunc(nops(gu2)/2))}):
      print(gu, `did not work out`):
    fi:
    fi:
   else

       if nops(gu2)>=2 then

        if  {seq(gu2[2*i1-1],i1=1..trunc(nops(gu2))/2)}<>{0} then
         print(gu, `did not work out`):
      fi:
      fi:
 fi:

   RETURN(gu):
  fi:
od:

FAIL:
end:


#Check34(m,NuP,Maxd): inputs an odd positive integer NuP (for the number of primes to check)
# to  test, and Maxd for the call of Check33 (make it big enough)
#and outputs pair [a_m,c_m], where a_m, c_m promised by Theorem 3.4 in the paper
#and checks that indeed the congruence, mod p^4 is right for all primes<=5 and less than MaxP. Try
#Check34(7,10,20);
Check34:=proc(m,NuP, Maxd) local am, cm, gu,X,k,i,Sm,k1,ku,p:

if not (type(m,integer) and m mod 2=1 and type(NuP,integer) and type(Maxd,integer) and NuP>0 and Maxd>0) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=Check33(k,1/2,3,m,Maxd) :

X:=gu[1]:

am:=gu[2][2]:
cm:=-4*subs(k=2,X):

for i from 3 to 3+NuP do
 p:=ithprime(i):

Sm:=add((-1)^k1*(4*k1+1)^m*(rf(1/2,k1)/k1!)^3,k1=0..(p-1)/2):

ku:=Sm-am*(p*(-1)^((p-1)/2) +p^3*En(p-3)) -p^3*cm:
ku:=ku mod p^4:
if ku<>0 then
 print([am,cm], `did not work out for the prime`, p):
 RETURN(FAIL):
fi:
od:

[am,cm]:
end: