######################################################################
##SAHI.txt: Save this file as  SAHI.txt                              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read SAHI.txt<Enter>                               #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Yukun Yao and Doron Zeilberger, Rutgers University      #
#[yao,zeilberg] at math dot rutgers dot edu                          #
######################################################################
 
#Created: 
 
print(`Created:  April, 2015. Renewed Aug. 2019`):
print(` This is SAHI.txt `):
print(`It is one of the packages that accompany the  article (in planning)  `):
print(`On Siddhartha Sahi's  Higher Correlation Inequalities`):
print(`by Yukun Yao and Doron Zeilberger`):

print(`and also available from Yao's and Zeilberger's websites`):
print(`It implements Siddhatha's Sahi's interesting article`):
print(`Higher Correlation Inequalities,  COMBINATORICA 28(2) (2008), 209-227, available from Sahi's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

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

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

ezraV:=proc()

if args=NULL then
 print(` The verbose procedures are:  RandI1v, RandIv(N,n,t,d)  `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  FixCar, IncFPS, IsPos1, NtoSet, randPP, RandI1, RandPos, RandPos1, SetToN `):
 print(` Zope `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: CheckC4,CheckC4m, CheckT1, CheckZ, Fplus,  MinT, RandI, Sid, Tay  `):
 print(` `):




elif nops([args])=1 and op(1,[args])=CheckC4 then
print(`CheckC4(N,n,d,[1/2,1/2,1/2],K): Checks Sahi's Conjecture 4 `):
print(`for random increasing functions (with coefficients from 0 to n),`):
print(`for the Boolean lattice with N elements`):
print(`degree d, and the specified m. It looks at the coefficients up to the K-th power. Try:`):
print(`CheckC4(3,2,3,[1/2,1/2,1/2],10);`):

elif nops([args])=1 and op(1,[args])=CheckC4m then
print(`CheckC4m(N,n,d,[1/2,1/2,1/2],K,L):  does CheckC4m(N,n,d,[1/2,1/2,1/2],K) L times. Try:`):
print(`CheckC4m(3,2,3,[1/2,1/2,1/2],10,100);`):

elif nops([args])=1 and op(1,[args])=CheckT1 then
print(`CheckT1(N,n,d,[1/2,1/2,1/2],K): Checks Sahi's Theorem 1 for random positive functions (with coefficients from 0 to n),`):
print(`for the Boolean lattice with N elements`):
print(`degree d, and the specified m. It looks at the coefficients up to the K-th power. Try:`):
print(`CheckT1(3,2,3,[1/2,1/2,1/2],10);`):

elif nops([args])=1 and op(1,[args])=CheckZ then
print(`CheckZ(d,t,N,a,K): Checks that if f(t) and g(t) are random polynomials of degree d with positive coefficients,`):
print(`then the  Maclaurin expansion of 1-(1-f(t))^a*(1-g(t))^(1-a) have non-negative coefficients up to order K. Try:`):
print(`CheckZ(3,t,30,1/2,30);`):

elif nops([args])=1 and op(1,[args])=FixCar then
print(`FixCar(L): inputs a list of length n+1, say, and outputs the numerical function that assigns to`):
print(`a subest of {1,..,n} L[|S|+1]. Try:`):
print(`FixCar([1,2]);`):

elif nops([args])=1 and op(1,[args])=Fplus then
print(`Fplus(F): given a function F, finds the function Fplus defined by equation (1) in Sahi's paper. Try`):
print(`Fplus([a,b]);`):

elif nops([args])=1 and op(1,[args])=IncFPS then
print(`IncFPS(L,t): given a list of lists each of length n+1, say, creates the`):
print(`function in the subsets of {1,...n}, such that each coeff. only depends on the cardinality `):
print(` where the coefficient of t^i goes by L[i] `):


elif nops([args])=1 and op(1,[args])=IsPos then
print(`IsPos(P,var): inputs a polynomial P in the list of variables var, and outputs true iff`):
print(`the coefficients are all non-negative. Try:`):
print(`IsPos(x^2+y^2,[x,y]);`):

elif nops([args])=1 and op(1,[args])=MinT then
print(`MinT(F,var,K): the set of coefficients, followed by the minimum, and 2nd-largest, of the Taylor expansion of the function F in the list of variables var`):
print(`up to order K. Try:`):
print(`MinT(1/(1-f1-f2),[f1,f2],10);`):

elif nops([args])=1 and op(1,[args])=NtoSet then
print(`NtoSet(N): the set corresponding to N, Try:`):
print(`NtoSet(3);`):


elif nops([args])=1 and op(1,[args])=randPP then
print(`randPP(t,d,N): a random polynomial in t (with constant term 0), of degree d with positive coefficients `):
print(`whose numerator and denominator are between 1 and N  Try:`):
print(`randPP(t,4,3); `):

elif nops([args])=1 and op(1,[args])=RandI then
print(`RandI(N,n,t,d): a random increasing polynomial-in-t-of-degree-d (that is 0 when t=0)`):
print(`function on the`):
print(` Boolean lattice where`):
print(` the difference between f(T)-f(S) is at most n when T covers S. `):
print(` Try: `):
print(` RandI(3,2,t,2); `):

elif nops([args])=1 and op(1,[args])=RandIv then
print(`RandIv(N,n,t,d): a verbose version of RandI(N,n,t,d): given as a list of pairs [Set,Value]. Only for humans. Try:`):
print(`RandIv(3,2,t,4);`):


elif nops([args])=1 and op(1,[args])=RandI1 then
print(`RandI1(N,n): a random increasing function on the Boolean lattice where`):
print(`the difference between f(T)-f(S) is at most n when T covers S.`):
print(`Try: `):
print(`RandI1(3,2);`):

elif nops([args])=1 and op(1,[args])=RandI1v then
print(`RandI1v(N,n): a verbose version of RandI1(N,n): given as a list of pairs [Set,Value]. Only for humans. Try:`):
print(`RandI1v(3,2);`):

elif nops([args])=1 and op(1,[args])=RandPos then
print(`RandPos(N,n,t,d): inputs a positive integer N and outputs a  random non-negative polynomial function in t, of degree d, drawn from {0,...,n}`):
print(`Try: RandPos(3,1,t,2);`):

elif nops([args])=1 and op(1,[args])=RandPos1 then
print(`RandPos1(N,n): inputs a positive integer N and outputs a  random non-negative function drawn from {0,...,n}`):
print(`Try: RandPos1(3,1);`):

elif nops([args])=1 and op(1,[args])=SetToN then
print(`SetToN(S): Given a set S of integers finds its code Try:`):
print(`SetToN({1,3});`):

elif nops([args])=1 and op(1,[args])=Sid then
print(`Sid(F,m): inputs a function F, on the subsets of {1,..,n},  given as a list of length 2^n, where`):
print(`the i-th item is the value of F at the subset of {1, ..,n} whose code is i, and a list of numbers between 0 and 1 (or left symbolic)`):
print(`of length n, outputs the expression given in Conjecture 4. Try:`):
print(`Sid([F00,F10,F01,F11],[m1,m2]);`):

elif nops([args])=1 and op(1,[args])=Tay then
print(`Tay(F,var,K): the Taylor expansion of the function F in the list of variables var`):
print(`up to order K. Try:`):
print(`Tay(1/(1-f1-f2),[f1,f2],10);`):

elif nops([args])=1 and op(1,[args])=Zope then
print(` Zope(F,G,m): 1-(1-F)^m*(1-G)^(1-m)`):
print(`Try:`):
print(`Zope(t,t^2,1/2); `):



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


#SetToN(S): Given a subset S of {1,...,N} finds its code between 1 and 2^N. Try:
#SetToN({1,3});
SetToN:=proc(S) local s:
1+add(2^(s-1), s in S):
end:

#NtoSet(N): the set corresponding to N, Try:
#NtoSet(3);
NtoSet:=proc(N) local r:
option remember:
if N=1 then
 RETURN({}):
fi:

r:=trunc(log[2](N-1)):

{r+1} union NtoSet(N-2^r):

end:

#Fplus(F): given a function F, finds the function Fplus defined by equation (1) in Sahi's paper. Try
#Fplus([a,b]);
Fplus:=proc(F) local N,S,i,gu,mu,mu1:
N:=log[2](nops(F)):

if not type(N, integer) then
 RETURN(FAIL):
fi:

gu:=[]:
for i from 1 to 2^N do
S:=NtoSet(i):
mu:=powerset(S):
gu:=[op(gu), add(F[SetToN(mu1)],mu1 in mu)]:
od:

gu:


end:


#Sid(F,m): inputs a function F, on the subsets of {1,..,n},  given as a list of length 2^n, where
#the i-th item is the value of F at the subset of {1, ..,n} whose code is i, and a list of numbers between 0 and 1 (or left symbolic)
#of length n, outputs the expression given in Conjecture 4. Try:
#Sid([F00,F10,F01,F11],[m1,m2]);
Sid:=proc(F,m) local i,s,S,n,gu,lu,j:

if not (type(m,list) and type(F,list)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

n:=nops(m):

if nops(F)<>2^n then
 print(`Bad input`):
 RETURN(FAIL):
fi:


S:=powerset({seq(i,i=1..n)}):

gu:=1:

for s in S do
lu:=1:

 for j from 1 to n do
   if member(j,s) then
     lu:=lu*m[j]:
   else
     lu:=lu*(1-m[j]):
   fi:
od:

gu:=gu*(1-F[SetToN(s)])^lu:

od:

1-gu:

end:


#FixCar(L): inputs a list of length n+1, say, and outputs the numerical function that assigns to
#a subest of {1,..,n} L[|S|+1]. Try:
#FixCar([1,2]);
FixCar:=proc(L) local n,i:
n:=nops(L)-1:
[seq(L[nops(NtoSet(i))+1],i=1..2^n)]:
end:


#IncFPS(L,t): given a list of lists each of length n+1, say, creates the
#function in the subsets of {1,...n}, such that each coeff. only depends on the cardinality
# where the coefficient of t^i goes by L[i]
IncFPS:=proc(L,t) local n,gu,j,i,lu:

n:=nops(L[1])-1:
gu:=[0$(2^n)]:

for i from 1 to nops(L) do
 lu:=FixCar(L[i]):
 gu:=[seq(gu[j]+lu[j]*t^i,j=1..nops(lu))]:
od:

gu:


end:


#RandPos1(N,n): inputs a positive integer N and outputs a  random non-negative function drawn from {0,...,n}
#Try: RandPos1(3,1);
RandPos1:=proc(N,n) local i,ra:
ra:=rand(0..n):
[seq(ra(),i=1..2^N)]:
end:

#RandPos(N,n,t,d): inputs a positive integer N and outputs a  random non-negative polynomial function in t, of degree d, drawn from {0,...,n}
#Try: RandPos(3,1,t,2);
RandPos:=proc(N,n,t,d) local i,ra,j:
ra:=rand(0..n):
[seq(add(ra()*t^j,j=1..d),i=1..2^N)]:
end:

#CheckT1(N,n,d,[1/2,1/2,1/2],K): Checks Sahi's Theorem 1 for random positive functions (with coefficients from 0 to n),
#for the Boolean lattice with N elements
#degree d, and the specified m. It looks at the coefficients up to the K-th power. Try:
#CheckT1(3,2,3,[1/2,1/2,1/2],10);
CheckT1:=proc(N,n,d,m,K) local t,gu,i:

gu:=RandPos(N,n,t,d):
gu:=Fplus(gu):
gu:=Sid(gu,m):

gu:=taylor(gu,t=0,K+1):

{seq(evalb(expand(coeff(gu,t,i))>=0),i=1..K)}:

end:


#RandI1(N,n): a random increasing function on the Boolean lattice where
#the difference between f(T)-f(S) is at most n when T covers S.
#Try:
#RandI1(3,2);
RandI1:=proc(N,n) local ra,f,i,j,gu,gu1,gu11,gu111,gadol,a:
ra:=rand(0..n):

f[{}]:=ra():

for i from 1 to N do
 gu:=choose({seq(j,j=1..N)},i):
  for gu1 in gu do
     gu11:={seq(gu1 minus {a}, a in gu1)}:
        gadol:=max(seq(f[gu111],gu111 in gu11)):
       f[gu1]:=gadol+ra():
   od:
od:

[seq(f[NtoSet(i)],i=1..2^N)]:
end:



#RandI(N,n,t,d): a random increasing polynomial-in-t-of-degree d function on the Boolean lattice where
#the difference between f(T)-f(S) is at most n when T covers S.
#Try:
#RandI(3,2,t,2);
RandI:=proc(N,n,t,d)  local i,i1,gu,gu1:

gu:=[0$(2^N)]:

for i from 1 to d do
 gu1:=RandI1(N,n):
 gu:=[seq(gu[i1]+gu1[i1]*t^i,i1=1..nops(gu))]:
od:

gu:

end:




#CheckC4(N,n,d,[1/2,1/2,1/2],K): Checks Sahi's Conjecture 4 
#for random increasing functions (with coefficients from 0 to n),
#for the Boolean lattice with N elements
#degree d, and the specified m. It looks at the coefficients up to the K-th power. Try:
#CheckC4(3,2,3,[1/2,1/2,1/2],10);
CheckC4:=proc(N,n,d,m,K) local t,gu,i:

gu:=RandI(N,n,t,d):
gu:=Sid(gu,m):

gu:=taylor(gu,t=0,K+1):

{seq(evalb(expand(coeff(gu,t,i))>=0),i=1..K)}:

end:




#RandI1v(N,n): a verbose version of RandI1(N,n): given as a list of pairs [Set,Value]. Only for humans. Try:
#RandI1v(3,2);
RandI1v:=proc(N,n) local gu,i:
gu:=RandI1(N,n):
[seq([NtoSet(i),gu[i]],i=1..nops(gu)) ]:
end:


#RandIv(N,n,t,d): a verbose version of RandI(N,n,t,d): given as a list of pairs [Set,Value]. Only for humans. Try:
#RandIv(3,2,t,2);
RandIv:=proc(N,n,t,d) local gu,i:
gu:=RandI(N,n,t,d):
[seq([NtoSet(i),gu[i]],i=1..nops(gu)) ]:
end:



#CheckC4m(N,n,d,[1/2,1/2,1/2],K,L): Checks Sahi's Conjecture 4 
#for random increasing functions (with coefficients from 0 to n),
#for the Boolean lattice with N elements
#degree d, and the specified m. It looks at the coefficients up to the K-th power. It does it for L cases
#Try:
#CheckC4m(3,2,3,[1/2,1/2,1/2],10,50);
CheckC4m:=proc(N,n,d,m,K,L) local gu,i:

for i from 1 to L do
 gu:=CheckC4(N,n,d,m,K):

  if gu<>{true} then
    RETURN(false):
  fi:
od:

true:

end:

#Zope(F,G,m): 1-(1-F)^m*(1-G)^(1-m)
Zope:=proc(F,G,m):
1-(1-F)^m*(1-G)^(1-m)
end:


#Tay(F,var,K): the Taylor expansion of the function F in the list of variables var
#up to order K. Try:
#Tay(1/(1-f1-f2),[f1,f2],10);
Tay:=proc(F,var,K) local gu,x,var1,lu,i,lu1:

if nops(var)=1 then
gu:=taylor(F,var[1],K+1):
gu:=add(normal(coeff(gu,var[1],i))*var[1]^i,i=0..K):
RETURN(gu):
fi:

var1:=[op(1..nops(var)-1,var)]:
x:=var[nops(var)]:

lu:=Tay(F,[x],K):
gu:=0:
for i from 0 to K do
lu1:=coeff(lu,x,i):
lu1:=Tay(lu1,var1,K):
gu:=expand(gu+lu1*x^i):
od:
gu:

end:



#MinT(F,var,K): the set of coefficients, followed by the minimum, and the second-largest,
#of the Taylor expansion of the function F in the list of variables var
#up to order K. Try:
#MinT(1/(1-f1-f2),[f1,f2],10);
MinT:=proc(F,var,K) local gu,gu1,i,j:
gu:={Tay(F,var,K)}:

for i from 1 to nops(var) do
gu:={seq(seq(coeff(gu1,var[i],j),j=0..K),gu1 in gu)}:
od:

gu,[min(gu), min(gu minus {0})]:

end:




#randPP(t,d,N): a random polynomial in t (with constant term 0), of degree d with positive coefficients 
#whose numerator and denominators are between 1 and N and denominator between 1 and M. Try:
#randPP(t,4,3);
randPP:=proc(t,d,N) local i:
add(rand(1..N)()/rand(1..N)()*t^i,i=1..d):
end:


#CheckZ(d,t,N,a,K): Checks that if f(t) and g(t) are random polynomials of degree d with positive coefficients,
#then the  Maclaurin expansion of 1-(1-f(t))^a*(1-g(t))^(1-a) have non-negative coefficients up to order K. Try:
#CheckZ(3,t,30,1/2,30);
CheckZ:=proc(d,t,N,a,K) local f,g,gu,i:
f:=randPP(t,d,N):
g:=randPP(t,d,N):
gu:=1-(1-f)^(1-a)*(1-g)^a:
gu:=taylor(gu,t=0,K+1):
if {seq(evalb(coeff(gu,t,i)>=0),i=1..K)}={true} then
 true:
else
false,[f,g]:
fi:

end:




#IsPos(P,var): inputs a polynomial P in the list of variables var, and outputs true iff
#the coefficients are all non-negative. Try:
#IsPos(x^2+y^2,[x,y]);
IsPos:=proc(P,var) local n,i,x,var1:

n:=nops(var):
x:=var[n]:


if n=1 then
for i from 0 to degree(P,x) do
 if coeff(P,x,i)<0 then
    RETURN(false):
 fi:
od:

RETURN(true):

fi:

var1:=[op(1..n-1,var)]:

for i from 0 to degree(P,x) do
 if not IsPos(coeff(P,x,i),var1) then
    RETURN(false):
 fi:
od:

true:

end: