######################################################################
##SahiProperty: Save this file as  SahiProperty.txt                  #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read SahiProperty.txt <Enter>                      #
##Then follow the instructions given there                           #
##Written by Yukun Yao and                                           # 
#Doron Zeilberger, Rutgers University ,                              #
#yao at math dot rutgers.edu DoronZeil at gmail dot com              #
######################################################################
 
#Created: 
 
print(`Created:  Sept. 2019`):
print(` This is SahiProperty.txt `):
print(`It is one of the packages that accompany the article (that may or may not come out) `):
print(` Epxerimenting with the Sahi Positivity property for certain formal power series`):
print(`by Yukun Yao 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: IsPosPol `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: CX, SahiE, SP `):
 print(` `):


elif nops([args])=1 and op(1,[args])=CX then
print(`CX(a,n): constructs a generic member of C[X} using a as in indexed variable. Try:`):
print(`CX(a,3);`):


elif nops([args])=1 and op(1,[args])=IsPosPol then
print(`IsPosPol(P,var): Given a polynomial P in a list of variables var, are the coefficients all non-negative?`):
print(`Try: `):
print(`IsPos(x+y,[x,y]); IsPos(x-y,[x,y]);`):


elif nops([args])=1 and op(1,[args])=SahiE then
print(`SahiE(L, var,alpha, K): The Sahi expression`):
print(`inputs a list L of linear expressions in a set of variables var, `):
print(`and a set of numbers or symbols alpha of the same lenght of L`):
print(`and a positive integer K outputs the first K coefficients of the Sahi-type expression`):
print(`1-(mul((1-L[i])^alpha[i],i=1..nops(L)) has a all positive coefficients checking up to degree K terms`):
print(` Try:  `):
print(`SahiE([x],[x],[1/2],10);`):
print(` SahiE([x,y],[x,y],[1/2,1/2],3);`):
print(` SahiE([x,y],[x,y],[a,b],3);`):
print(`SahiE([x,x+y],[x,y],[1/2,1/2],10);`):
print(`For the first 4 terms for Sahi theorm for B_2`):
print(`SahiE([a,a+b,a+c,a+b+c+d],[a,b,c,d],[1/4,1/4,1/4,1/4],4);`):
print(`For the first 4 terms of Lemma 16 in Sahi's paper:`):
print(`SahiE([a,a+b+c,a+c,a+b+c+d+e],[a,b,c,d,e],[1/4,1/4,1/4,1/4,1/4],4);`):

elif nops([args])=1 and op(1,[args])=SP then
print(`SP(L, var,alpha, K): Testing the Sahi property `):
print(`inputs a list L of linear expressions in a set of variables a, and a positive integer K, checks`):
print(`whether 1-(mul(1-L[i],i=1..nops(L))^alpha has a all positive coefficients checking up to degree K terms`):
print(`Try:  `):
print(`SP([x],[x],1/2,10);`):
print(` SP([x,y],[x,y],1/2,10);`):
print(`SP([x,x+y],[x,y],1/2,10);`):
print(`For the Sahi theorm for B_2`):
print(`SP([a,a+b,a+c,a+b+c+d],[a,b,c,d],1/4,10);`):
print(`(note that 1/4 seems to be best possible)`):
print(`For and empirical check up to 20 terms of Lemma 16 in Sahi's paper:`):
print(`SP([a,a+b+c,a+c,a+b+c+d+e],[a,b,c,d,e],1/4,10);`):
print(`(note that 1/4 is best possible)`):

print(``):


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




#IsPosPol(P,var): Given a polynomial P in a list of variables var, are the coefficients all non-neatgive?
#Try:
#IsPos(x+y,[x,y]); IsPos(x-y,[x,y]);
IsPosPol:=proc(P,var) local P1,gu,i,i1:

if not type(P,`+`) then
 gu:= P/mul(var[i]^degree(P,var[i]),i=1..nops(var)):
 if not type(gu, numeric) then
     RETURN(FAIL):
 fi:
 if gu<0 then
   RETURN(false):
 fi:
fi:

for i from 1 to nops(P) do
P1:=op(i,P):

 gu:=normal(P1/mul(var[i1]^degree(P1,var[i1]),i1=1..nops(var))):

 if not type(gu, numeric) then
     RETURN(FAIL):
 fi:
 if gu<0 then
   RETURN(false):
 fi:

od:

true:
end:

#SP(L, var,alpha, K): Testing the Sahi property 
#inputs a list L of linear expressions in a set of variables a, and a posotive integer K, checeks
#whether 1-(mul(1-L[i],i=1..nops(L))^alpha has a all positive coefficients checking up to degree K terms
#Try: 
#SP([[x],[x],1/2,10);
# SP([[x,y],[x,y],1/2,10);
#SP([[x,x+y],[x,y],1/2,10);
#For the Sahi theorm for B_2
#SP([a,a+b,a+c,a+b+c+d],[a,b,c,d],1/4,10);
#(note that 1/4 is best possible)
#For and empirical check up to 20 terms of Lemma 16 in Sahi's paper:
#SP([a,a+b+c,a+c,a+b+c+d+e],[a,b,c,d,e],1/4,10);
#(note that 1/4 is best possible)

SP:=proc(L, var,alpha, K) local gu,i,t,gu1:
gu:=1-mul(1-L[i],i=1..nops(L))^alpha:
gu:=subs({seq(var[i]=t*var[i],i=1..nops(var))},gu):

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

for i from 1 to K do
 gu1:=expand(coeff(gu,t,i)):
 if not IsPosPol(gu1,var) then
   print(`The Sahi property fails at degree`, i):
  RETURN(false):
 fi:
od:
true:

end:



#CX(a,n): constructs a generic member of C[X} using a as in indexed variable. Try:
#CX(a,3);
CX:=proc(a,n) local mu,mu1,mu11,i,gu:
mu:=convert(powerset(n),list):


gu:=[]:

for i from 1 to nops(mu) do
 mu1:=powerset(mu[i]):

 gu:=[op(gu),add(a[mu11],mu11 in mu1)]:
od:
 

gu,[seq(a[mu[i]],i=1..nops(mu))]:

end:



#SahiE(L, var,alpha, K): The Sahi expression
#inputs a list L of linear expressions in a set of variables var, 
#and a set of numbers or symbols alpha of the same lenght of L
#and a positive integer K outputs the first K coefficients of the Sahi-type expression
#1-(mul((1-L[i])^alpha[i],i=1..nops(L)) has a all positive coefficients checking up to degree K terms
#Try: 
#SahiE([[x],[x],[1/2],10);
# SahiE([[x,y],[x,y],[1/2,1/2],3);
# SahiE([[x,y],[x,y],[a,b],3);
#SahiE([[x,x+y],[x,y],[1/2,1/2],10);
#For the Sahi theorm for B_2
#SahiE([a,a+b,a+c,a+b+c+d],[a,b,c,d],[1/4,1/4,1/4,1/4],3);
#For and empirical check up to 20 terms of Lemma 16 in Sahi's paper:
#SahiE([a,a+b+c,a+c,a+b+c+d+e],[a,b,c,d,e],[1/4,1/4,1/4,1/4,1,4],3);

SahiE:=proc(L, var,alpha, K) local gu,i,t,gu1:
if nops(L)<>nops(alpha) then
print(L, `and `, alpha, `should be the same length `):
 RETURN(FAIL):
fi:
gu:=1-mul((1-L[i])^alpha[i],i=1..nops(L)):
gu:=subs({seq(var[i]=t*var[i],i=1..nops(var))},gu):

gu:=taylor(gu,t=0,K+1):
expand([seq(coeff(gu,t,i),i=1..K)]):

end: