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


print(`Written : Dec. 2021- Jan. 2022  `):
print():
print(`This is CritVal.txt, A Maple package`):
print(`accompanying Tewodros Amdeberhan and Doron Zeilberger's  article: `):
print(` An elegant Multi-Integral that implies an even more elegant determinant identity of Dougherty and McCammond  `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/CritVal.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
 



ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(` F, L, R `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(`  CritVal.txt: A Maple package for  confirming the statements in the Amdeberhan-Zeilberger paper`):

print(`The MAIN procedures are`):
print(` CheckCauchy,  CheckEq1, CheckEq2, CheckEq3L, CheckEq3R, CheckEq4, CheckEq5, CheckEq6 `):



elif nargs=1 and args[1]=CheckCauchy then
print(`CheckCauchy(n): checks Cauchy's alternant formula for n variables. Try:`):
print(`CheckCauchy(3);`):

elif nargs=1 and args[1]=CheckEq4 then
print(`CheckEq4(a): checks Eq. (4)  for the list of integers a. Try:`):
print(`CheckEq4([2,3,4]);`):


elif nargs=1 and args[1]=CheckEq5 then
print(`CheckEq5(a): Checks Eq (5) for the list of integers a (and symbolic z). Try:`):
print(`CheckEq5([1,2,3]);`):

elif nargs=1 and args[1]=CheckSomething then
print(`CheckSomething(a,i0): Checks something obvious permutation pi of {1,...,n-1}. Try:`):
print(`CheckSomething([2,3,4],1);`):


elif nargs=1 and args[1]=CheckEq6 then
print(`CheckEq6(n,A,z): Checks Eq. (6) for dimension n. It reuturns true and the actual symbolic expression of the signed sums of regions.Try:`):
print(`CheckEq6(4,A,z);`):

elif nargs=1 and args[1]=CheckEq1 then
print(`CheckEq1(a,b): Checks the identity in Theorem A (Eq. 1) for the list of integers a, and the positive integer b.Try:`):
print(`CheckEq1([2,3,4],3);`):


elif nargs=1 and args[1]=CheckEq2 then
print(`CheckEq2(a): Checks the identity in (Eq2), i.e. Theorem B for the list of integers a. Try:`):
print(`CheckEq2([2,3,4]);`):


elif nargs=1 and args[1]=CheckEq3L then
print(`CheckEq3L(n): Checks that the recurrence (Eq3) is satisfied by the integrand of the left side of Eq.(1). Try:`):
print(`Try: CheckEq3L(3);`):


elif nargs=1 and args[1]=CheckEq3R then
print(`CheckEq3R(n): Checks that the recurrence (Eq3) is satisfied by the right side of Eq.(1). Try:`):
print(`Try: CheckEq3R(3);`):

elif nargs=1 and args[1]=F then
print(`F(z,x,a,b): Given a list,a, of non.-neg. pos. integers finds the integrand of the left side of the identity. Try:`):
print(`F([z1,z2],[x1,x2],[a1,a2],b);`):

elif nargs=1 and args[1]=L then
print(`L(a,z,b): The left side of the identity in (Eq1)  in Theorem A. Try:`):
print(`L([2,3,4],[z1,z2,z3],3);`):


elif nargs=1 and args[1]=R then
print(`R(a,z,b): The right side of the identity in (Eq1) in Theorem A. Try:`):
print(`R([2,3,4],[z1,z2,z3],3);`):

else
 print(`There is no such thing as`, args):

fi:


end:

with(linalg):
with(combinat):


#V(x): The Vandermone of the list x. Try; V([a,b,c]);
V:=proc(x) local i,j: mul(mul(x[j]-x[i],j=i+1..nops(x)),i=1..nops(x)):end:
	   

#F(z,x,a,b): Given a list,a, of non.-neg. pos. integers finds the integrand of the left side of the identity. Try:
#F([z1,z2],[x1,x2],[a1,a2],b);
F:=proc(z,x,a,b) local n,i,j,k:
n:=nops(a):
mul(x[i]^b,i=1..n)*mul( (mul(x[j]-z[k],j=1..n))^a[k],k=1..n):
end:


#CheckEq4(a): Checks Eq4 in the paper
CheckEq4:=proc(a) local n,z,x,i,j,L,gu:
n:=nops(a):
L:=F([seq(z[i],i=1..n)], [seq(x[i],i=1..n)],a,0):



gu:=0:
for i from 1 to n do
gu:=gu+(-1)^i*diff( L*mul(x[i]-z[j],j=1..n),x[i])*V([seq(x[j],j=1..i-1), seq(x[j],j=i+1..n)]):
od:
evalb(normal((-1)^n*gu/L/V([seq(x[i],i=1..n)]))=add(a[i],i=1..n)+n):
end:



#L(a,z,b): The left side of the identity in Eq1. Try:
#L([2,3,4],[z1,z2,z3],3);
L:=proc(a,z,b) local x,n,i,gu:
n:=nops(a):

gu:=F([seq(z[i],i=1..n)], [seq(x[i],i=1..n)],a,b)*V([seq(x[i],i=1..n)]):

for i from 1 to n do
gu:=int(gu,x[i]=0..z[i]):
od:
factor(gu):
end:


#R(a,z,b): The right side of the identity in Eq1
R:=proc(a,z,b) local n,i,j:
n:=nops(a):
(-1)^add(i*a[i],i=1..n)*
mul(mul((z[j]-z[i])^(a[j]+a[i]+1),j=i+1..n),i=1..n)*mul(z[j]^(a[j]+b+1),j=1..n)*b!*mul(a[j]!,j=1..n)/(add(a[j],j=1..n)+n+b)!:
end:


#CheckEq1(a,b): Checks the identity in Eq1 for the list of integers a,  and the positive integer b.Try:
#CheckEq1([2,3,4],3);
CheckEq1:=proc(a,b) local n,i,z :
n:=nops(a):

evalb(normal(L(a,[seq(z[i],i=1..n)],b)/R(a,[seq(z[i],i=1..n)],b))=1):
end:



#CheckEq5(a): Check Eq. (5) Try:
#CheckEq5([1,2,3]);
CheckEq5:=proc(a) local z,n,y,i,j,L,gu,gu1,lu:
n:=nops(a):
L:= mul(mul(y[i]-z[j],i=1..n-1)^a[j],j=1..n)*V([seq(y[i],i=1..n-1)]):



gu:=[]:

for i from 1 to n do
gu1:=L:

for j from 1 to i-1 do
 gu1:=int(gu1,y[j]=0..z[j]):
od:

for j from i+1 to n do
 gu1:=int(gu1,y[j-1]=0..z[j]):
od:
gu1:=factor(gu1):
gu:=[op(gu),(-1)^(i-1)*gu1]:

od:

gu:

lu:=L:

for i from 1 to n-1 do
lu:=int(lu,y[i]=z[1]..z[i+1]):
od:
lu:=factor(lu):

evalb(normal(factor(convert(gu,`+`))/lu)=1):

end:





#pza(Z,z,a): Integral of (w-z[i])^a[1]*...*(w-z[n])^a[n] w.r.t. w, from 0 to Z
pza:=proc(Z,z,a) local i,w,n:

n:=nops(a):
expand(int(mul((w-z[i])^a[i],i=1..n),w=0..Z)):
end:

#T3(z,n): The expression in Theorem 1
T3:=proc(z,a) local i,j,gu,Z,n:
n:=nops(a):
gu:=pza(Z,z,a):

factor(det([seq(   [seq(diff(subs(Z=z[j],gu),z[i]),j=1..n)],i=1..n)])):
end:

R3:=proc(z,a) local n,i,j,k: 
n:=nops(a):
factor(mul(a[i]!,i=1..n)/add(a[i],i=1..n)!*mul((-z[j])^a[j],j=1..n)*mul((mul(z[k]-z[j],k=1..j-1)*mul(z[k]-z[j],k=j+1..n))^a[j],j=1..n)):
end:


CheckEq2:=proc(a) local z:
evalb(normal(T3(z,a)/R3(z,a))=1) :
end:

#CheckSomething(a,i0): Checks some property
#CheckSomething([2,3,4],[1,2]);
CheckSomething:=proc(a,i0) local F,y,z,n,L1,L2,P,GU,i,k,pi,j:
n:=nops(a):
F:=mul(mul(y[j]-z[i],j=1..n-1)^a[i],i=1..n)*V([seq(y[i],i=1..n-1)]):

L1:=F:

for i from 1 to i0-1 do
L1:=normal(int(L1,y[i]=0..z[i])):
od:
for i from i0 to  n-1 do
L1:=normal(int(L1,y[i]=0..z[i+1])):
od:

L1:=factor(L1):
GU:=[]:

P:=permute([seq(j,j=1..i0-1),seq(j,j=i0+1..n)]):

for k from 1 to nops(P) do
pi:=P[k]:
L2:=F:
for i from 1 to nops(pi) do
L2:=normal(int(L2,y[i]=0..z[pi[i]])):
od:
L2:=factor(L2):
GU:=[op(GU),[pi,normal(L2/L1)]]:
od:
GU:
end:







#sgn(pi): The sign of the permutation pi
sgn:=proc(pi) local i,j,c:
c:=1:

for i from 1 to nops(pi) do
 for j from i+1 to nops(pi) do
 if pi[i]>pi[j] then
  c:=-c:
 fi:
 od:
od:
c:
end:


#CheckEq6: Checks Eq. (6) for n.Try:
#CheckEq6(4);
CheckEq6:=proc(n,A,z) local F,i,gu,M1,pi,M2,i1,j,G,M3:
gu:=permute(n-1):

M1:=add((-1)^(i-1)*mul(A[i1][0,z[i1]],i1=1..i-1)*mul(A[i1][0,z[i1+1]],i1=i..n-1),i=1..n):



F:=0:

for pi in {op(gu)} do

M2:=subs({seq(A[i1]=A[pi[i1]],i1=1..n-1)},M1):
F:=F+sgn(pi)*M2:
od:


for i from 1 to n-1 do
 for j from 1 to n do
  F:=subs(A[i][0,z[j]]=A[i][0,z[1]]+A[i][z[1],z[j]],F):
 od:
od:

F:=expand(F):

for i from 1 to n-1 do
F:=subs(A[i][z[1],z[1]]=0,F):
od:


M2:=mul(A[i][z[1],z[i+1]],i=1..n-1):
G:=0:

for pi in {op(gu)} do

M3:=subs({seq(A[i1]=A[pi[i1]],i1=1..n-1)},M2):
G:=G+sgn(pi)*M3:
od:


[evalb(F-G=0),F]:



end:


CheckEq3L:=proc(n) local i,j,x,z,smol,yemin:
smol:=mul(x[i],i=1..n):

yemin:=0:

for i from 1 to n do
yemin:=yemin+mul(z[j],j=1..n)/z[i]*mul(x[j]-z[i],j=1..n)/(mul(z[j]-z[i],j=1..i-1)*mul(z[j]-z[i],j=i+1..n)):
od:

yemin:=yemin+mul(z[j],j=1..n):
yemin:=normal(yemin):
evalb(smol=yemin):
end:


CheckEq3R:=proc(n) local i,j,a,z,smol,yemin,b,i1:

smol:=(b+1)/(n+b+1+ add(a[i],i=1..n))*mul(z[j],j=1..n):

yemin:=0:

for i from 1 to n do
yemin:=yemin-z[i]*(a[i]+1)/(n+b+1+ add(a[i1],i1=1..n))*mul(z[j],j=1..n)/z[i]:
od:

yemin:=yemin+mul(z[j],j=1..n):
yemin:=normal(yemin):
evalb(smol=yemin):
end:

CheckCauchy:=proc(n) local x,z,i,j,smol,yemin:
smol:=factor(det([seq([seq(1/(x[j]-z[i]),j=1..n)],i=1..n)])):
yemin:=factor(mul(mul((z[i]-z[j])*(x[j]-x[i]),j=i+1..n),i=1..n)/mul(mul(x[j]-z[i],i=1..n),j=1..n)):
evalb(normal(smol/yemin)=1):
end: