#file MultInt.txt

#last update: Wed, Oct 22, 2014 11:38 PM.

# This is  MultInt: a Maple package for MULTiple INTegration
# of proper-hyperexponential functions by the continuous version of
# the multi-WZ method.

# AUTHOR: Akalu Tefera
#        Grand Valley State University, Math Dept.
#        TeferaA@gvsu.edu

# MultInt is the updated and generalized version of the package
# TRIPLE_INTEGRAL which was written by Doron Zeilberger and MultInt is the
# improved version of the packages, Mint3 and Mint5.

# CAUTION: this version of MultInt is for Maple 7 and above.

# If you have any questions or want to report some bugs please contact the author.

# ACKNOWLEDGMENT: Many thanks to Doron Zeilberger, my thesis advisor,  for his
# helpful suggestions and valuable support. 
#         

# ------------------------------------------------------------------------------  
# `MultInt/freeintgrand`(intgrand, vars) return true if integrand does not contain all
# the variables in vars. A message is printed in this case.

`MultInt/freeintgrand`:=proc(intgrand1,vars)

local i,intgrand:

intgrand:=normal(simplify(intgrand1)):

if not type(vars,list) then ERROR(`expecting list`) fi:

 if member(false,convert([seq(has(intgrand,vars[i]),i=1..nops(vars))],set)) then

   RETURN(true)

 fi:

 RETURN(false):

end: 


# ------------------------------------------------------------------------------  
# `MultInt/checkpoly` : checks if input is list of poly of length len.
#            returns list of polys by replacing 0 polys by 1 polys
#            if the list contain non-poly then displays ERROR msg
  

`MultInt/checkpoly`:=proc(poly,len)

if nops(poly)= len then

 RETURN(subs(0=1,poly))

fi:

ERROR(cat(convert(poly,string),` must contain `, 

          convert(len,string),` polynomials.`)):

end:


#---------------------------------------------

`MultInt/extrctauxvars`:=proc(f,intnvars)

local vars:

vars:=indets(f):

vars:=select(type,vars,name):

vars:=vars minus intnvars:

vars:=convert(vars,list):

vars:

end:

#------------------------------------------------------------------------------
# `MultInt/notdistvars` returns true if the variable list contains multiple elements, 
# and false otherwise.
 
`MultInt/notdistvars`:=proc(vars)

if nops(convert(vars,set)) = nops(vars) then

 RETURN(false):

else

 RETURN(true):

fi:

end: 

#----------------------------------------------------------------------
# checks the equality of the sizes of the lists lst1 and lst2

`MultInt/equalsize`:=proc(lst1,lst2)

if nops(lst1) <> nops(lst2) then 

 RETURN(false)

else
  
 RETURN(true)

fi:

end:



#------------------------------------------------------------------------------
`MultInt/structurelst`:=proc(ansatz1,shvars)

local ansatz, strlst, mono, i, j, k, vfactrs, factrs, indx, v:

ansatz:=convert({op(ansatz1)},list):

strlst:=[]:

for i from 1 to nops(ansatz) do

 mono:=ansatz[i]:
 
 factrs:=factors(mono):
 
 if factrs[1]<>1 then
 
   ERROR(cat(` Invalid list `, convert(ansatz,string))):
   
 fi:
 
 factrs:=op(2,factrs):

 indx:=[seq(0,k=1..nops(shvars))]:

 for j from 1 to nops(shvars) while nops(factrs)<>0 do
  
  v:=shvars[j]:
   
  vfactrs:=select(has,factrs,v):
  
  if nops(vfactrs)=1 then
  
   vfactrs:=op(vfactrs):
   
   indx:=subsop(j=op(2,vfactrs),indx):
   
   factrs:={op(factrs)} minus {vfactrs}:
   
   factrs:=convert(factrs,list):
      
  elif nops(vfactrs)>1 then
  
   ERROR(cat(`invalid list `, convert(ansatz,string))):
   
  fi:
 
 od:

 strlst:=[op(strlst),indx]:
     
od:

strlst:={op(strlst)} minus {[seq(0,i=1..nops(shvars))]}:
  
strlst:=convert(strlst,list):

strlst:


end: 

#--------------------------------------------------------------------------
`MultInt/ver3or2`:=proc()

local st,ver,pos,subst,verset:

st:=interface(version):

ver:={}:

pos:=searchtext(`,`,st):

while pos<>0 do
  
 subst:=substring(st,1..pos-1):

 ver:=ver union {subst}:

 st:=substring(st,pos+1..length(st)):

 pos:=searchtext(`,`,st):

 if pos=0 then
  
  ver:=ver union {st}

 fi:

od:

verset:={` Maple V`,`Maple V`,` Release 3`, `Release 3`,`Release 2`,

        ` Release 2`,`TTY Iris`,` TTY Iris`}:
        
if nops(verset intersect ver) >=2 then

 true

else

 false

fi:


end:


#----------------------------------------------------------------------------
# `MultInt/makindx` computes, given a list of nonnegative integers, a list of indices 
# s.t. for each entry index[i] <= min(degs[i],
#                  tdeg-index[nops(index)]-sum(indx[j],j=1..(i-1))) holds 
# (for all i <= nops(degs)).

`MultInt/makindx`:=proc(degs,tdeg)

local indx,indxset,i,j,t:

t:=[list(nonnegint),nonnegint]:

if not type([args],t)  then

  ERROR(cat(`expecting an arg of type`,

  `list of integers(>=0), integer(>=0)`))

fi:

if member(true,convert( map(proc(x) x<0 end,[op(degs),tdeg]),set)) then

   ERROR(cat(`expecting an arg of type`,
   
   `list of integers(>=0), integer(>=0)`))

fi:

indx := [seq(0,i=1..nops(degs))]:

indxset := {}:

while indx[nops(indx)] <= degs[nops(degs)] do

 indxset:=indxset union {indx}:

 # increment index      

 i := 1:

  while indx[i]>= min(degs[i],tdeg-indx[nops(indx)]-sum('indx[j]','j'=1..(i-1))) 

         and (i<nops(degs)) do

    indx:=subsop(i=0,indx):

    i:=i+1:

  od:  

  indx:=subsop(i=indx[i] + 1,indx):

od:  

convert(indxset,list):

end:    

#-------------------------------------------------------------------------------
`MultInt/makindx1`:=proc(n)

local lst,i,t:

t:=[posint]:

if not type([args],t) then 

 ERROR(` expects an argument of type positive integer`):

fi:

lst:=[]:

for i from 1 to n-1 do

 lst:=[op(lst),[i,i+1]]:

od:

lst:

end:

# ---------------------------------------------------------------------------
# `MultInt/makindx2` computes, given a list of nonnegative integers, a list
#   of indices s.t. for each entry index[i] <= orders[i] holds 
#   (for all i <= nops(orders).
 
`MultInt/makindx2`:=proc(orders)

local indx,indxset,i,t:

t:=[list(nonnegint)]:

if not type([args],t)  then

  ERROR(` expecting an arg of type: list of integers(>=0)`)

fi:

if nops(args)=0 then ERROR(` invalid input `) fi:

if member(true,convert(map(proc(x) is(x<0) end,args),set)) then

   ERROR(`expecting an arg of type:list of positive integers`)

fi:

indx := [seq(0,i=1..nops(orders))]:

indxset := {}:

while indx[nops(indx)] <= orders[nops(orders)] do

 indxset:=indxset union {indx}:

 # increment index      

 i := 1:

 while indx[i]=orders[i] and  i<nops(orders) do

  indx:=subsop(i=0,indx):
  
  i:=i+1:

 od:  

 indx:=subsop(i=indx[i] + 1,indx):
   
od:  

indxset:= indxset minus {[seq(0,i=1..nops(orders))]}:

convert(indxset,list):

end: 
  

#------------------------------------------------------------------------------
`MultInt/chkD`:=proc(expr)

local result:

 if expr=1 then  
 
  result:=true 
 
 elif type(expr,`^`) and abs(op(1,expr)) = 1 and (op(2,expr) mod 2 ) = 0 then
 
   result:=true:
   
 else
 
  result:=false:
  
 fi:
 
result:
  
end:

#----------------------------------------------------------------------------
`MultInt/sym`:=proc(f1,intnvar)

local t,indxlst,i,indx,siderel,D1,f, result:

if nargs > 2 then
 
ERROR(` wrong number of arguments` );

fi:

t:=[algebraic, list(name)]:

f:=factor(simplify(f1)):

if not type([args],t) then

 ERROR(cat(` Expecting argument of type:\n`,

 `algebraic expression, list of var. names`)):

fi:

indxlst := `MultInt/makindx1`(nops(intnvar)):

result:=true:

for i from 1 to nops(indxlst) while result=true do

 indx:=indxlst[i]:

 siderel := {intnvar[indx[1]]=intnvar[indx[2]],intnvar[indx[2]]=intnvar[indx[1]]}:

 D1:=simplify(expand(eval(subs(siderel,f))/f)):
 
 result:=`MultInt/chkD`(D1):
 
od:

result:

end:

MultInt[sym]:=proc()

`MultInt/sym`(args)

end:

#--------------------------------------------------------------------------
`MultInt/esp`:=proc(n,vars)

local i,j,indxlst,s,indx,t:

t:=[posint,list(name)]:

if nargs > 2 then 

 ERROR(` wrong number of arguments`):
 
fi:

if not type([args],t) then
 
  ERROR(cat( `Expecting an argument of type:\n`,

  ` positive integer, variable names`)):
  
fi:
  
if nops(vars) < n then RETURN(0) fi:

with(combinat):

indxlst:=choose(vars,n):

s:=0:

for i from 1 to nops(indxlst) do
 
 indx:=indxlst[i]:
 
 s:= s + product(indx[j],j=1..nops(indx)):
 
od:

s:

end:

MultInt[esp]:=proc()

`MultInt/esp`(args):

end:


#------------------------------------------------------------------------------
# given a integrand f, with integration variables intnvars extracts the
# denominator of  lcm(D_xj (f)/f) 
 
`MultInt/extrctdenom`:=proc(f,intnvars)

local result,denset,i,den:

denset:={}:

for i from 1 to nops(intnvars) do

 result:=diff(simplify(log(f)),intnvars[i]):  

 result:=normal(result):

 den:=denom(result):

 denset:=denset union {den}:

od:

den := lcm(op(denset)):

den:= factor(den):

den:

end:


#-------------------------------------------------------------------------------
# Given f1, `MultInt/extrctP` extracts  Polynomial expression in the variables:
# mainvars and intnvars from f1:

`MultInt/extrctP`:=proc(f1,mainvars,intnvars,axvars)

local p, f, expr, i, axvars1:

f:=normal(convert(f1,factorial)):

p:=factor(numer(f)):  f:=1/denom(f):

if not type(p,polynom(anything,[op(intnvars),op(mainvars)])) then

 axvars1:={op(axvars)} minus {op(mainvars)}:

 for i from 1 to nops(axvars1) do

   expr:=select(has,``*p,axvars1[i]):
   
   if expr<>NULL then
   
     p:=factor(simplify(p/expr)):
 
     f:=factor(f*expr):
    
   fi:

 od:

 if not type(p,polynom(anything,[op(intnvars),op(mainvars)])) then
  
   expr:=select(type,``*p,polynom(anything,[op(intnvars),op(mainvars)])):
        
   if expr<>NULL and has(expr,[op(intnvars),op(mainvars)])  then
    
    f:=factor(f*factor(simplify(p/expr))):

    p:=expr:
    
   elif expr<>NULL then
   
    p:=factor(simplify(p/expr)):
 
    f:=factor(f*expr):
    
   fi:

  fi:

fi:

if not type(p,polynom(anything,[op(intnvars),op(mainvars)])) then

  f:=factor(f*p):
    
  p:=1:
  
fi:
 
[p,f]:

end:

#------------------------------------------------------------------------------
# `MultInt/fa1` and `MultInt/fa2`
# - find the functional equation for the rational functions
# R_i  satisfying, for a given proper-hyperexponential function f, 
# the functional equation: Operator(f) = D_y1 (R_1 f) +  D_y2 (R_2 f) + ...

`MultInt/fa1`:=proc(p,f, mainvar,intnvars,mainorder)

local i,gu:

gu:=p:

for i from 1 to mainorder do

 gu:=gu+(cat(_a_,i))*subs(mainvar=mainvar+i,p)*simplify(subs(mainvar=mainvar+i,f)/f):

od:

gu:=gu - sum('cat(YR_,i)','i'=1 .. nops(intnvars)) 

       - sum('(cat(_R_,i)) * normal(diff(simplify(log(f)),intnvars[i]))',
   
            'i'=1..nops(intnvars)):

gu:

end:

#------------------------------------------------------------------------------

`MultInt/fa2`:=proc(p,f, mainvars,intnvars,mainords_or_strlst)

local i,j,k,gu,indxlst,indx,cfindx:

gu:=p:

if type(mainords_or_strlst,list(nonnegint))  and nops(mainords_or_strlst)>0

   then

 indxlst:=`MultInt/makindx2`(mainords_or_strlst):
 
else

 indxlst:=mainords_or_strlst:
 
fi:

for i from 1 to nops(indxlst) do

 indx:= op(i,indxlst):
 
 cfindx:=(cat(seq(op(j,indx),j=1..nops(indx)))): 
 
 gu:=gu +
 
 (cat(_a_,cfindx))*subs({seq(mainvars[j]=mainvars[j]+op(j,indx),j=1..nops(indx))},p)*
     
  simplify(subs({seq(mainvars[k]=mainvars[k]+op(k,indx),k=1..nops(indx))},f)/f):

od:

gu:=gu - sum('cat(YR_,i)','i'=1 .. nops(intnvars)) 

      - sum('(cat(_R_,i)) * normal(diff(simplify(log(f)),intnvars[i]))',
   
            'i'=1..nops(intnvars)):

gu:

end:

#------------------------------------------------------------------------------
# `MultInt/pashet1` and `MultInt/pashe2`: writes a given operator in normal form

`MultInt/pashet1`:=proc(p,N,mainvar,intnvars,fnam)

local i,gu,p1,deg:

p1:=collect(p,N):

gu:=0:  deg:=degree(p,N):

if p1=0 and member(deg,{infinity,-infinity}) then deg:=0: fi:

for i from 0 to  deg do
   
 gu:=gu+factor(coeff(p1,N,i))*fnam(mainvar + i,op(intnvars)):

od:

RETURN(gu):

end:

#----------------------------------------------------------------------------
`MultInt/pashet2`:=proc(p,mainvars,intnvars,fnam)

local i,j,factrs,indx,gu,pcffs,v,t,pos:

v:=[seq(cat(_N_,i),i=1..nops(mainvars))]:

pcffs:=[coeffs(collect(p,v,'distributed'),v,'t')]:

t:=[t]:

gu:=0: 

for i from 1 to nops(pcffs) do
 
 factrs:=op(2,factors(t[i])):
 
 if nops(factrs)<>0 then
  
  indx:=[seq(mainvars[j],j=1..nops(mainvars))]:
  
  for j from 1 to nops(factrs) do
  
   pos:='pos':
    
   if member(op(1,op(j,factrs)),v,'pos') then
     
    indx:=subsop(pos=op(2,op(j,factrs)) + op(pos,mainvars),indx):
    
   fi:
  
  od:
    
  gu:=gu+factor(op(i,pcffs))*fnam(op(indx),op(intnvars)):
     
 else
 
  gu:= gu + factor(op(i,pcffs))*fnam(op(mainvars),op(intnvars)):
  
fi:

od:

RETURN(gu):

end:

#------------------------------------------------------------------------------
# `MultInt/cv`:  given a functional equation, eq, in 
#             Ri,YRi(=D_yi(Ri)),i=1..nops(intnvars), and some proposed 
#             Change of dependent Variables(and hence called cv)
#       OldRi:=NewRi*Si, i=1..nops(intnvars).
#
#      The output is the functional equation in the NewRi,i=1...,nops(intnvars). 


`MultInt/cv`:=proc(S,eq,intnvars)

local gu,i:

gu:=eq:

for i from 1 to nops(S) do

 gu:=subs( (cat(_R_,i)) = (cat(_R_,i))*S[i],gu):

 gu:=subs((cat(YR_,i)) = S[i]*(cat(YR_,i))+diff(S[i],intnvars[i])*(cat(_R_,i)),gu):

od:

 `MultInt/nor`(gu):

end:


#------------------------------------------------------------------------------
# `MultInt/nor`: normalizes each operands of bitu

`MultInt/nor`:=proc(bitu)

local i,gu:
 
gu:=0:

if type(bitu,`+`) then

 for i from 1 to nops(bitu) do 

  gu:=gu+normal(op(i,bitu)):

 od:

else

 gu:=normal(bitu):

fi:

 gu:

end:

#-----------------------------------------------------------------------------

`MultInt/finddeg`:=proc(eq,Ri,YRi,x)

local deg1,deg2,deg3,degset:

degset:={infinity,-infinity}:

deg1:=degree(collect(eq,x),x):

deg2:=degree(collect(coeff(collect(eq,Ri),Ri),x),x):

deg3:=degree(collect(coeff(collect(eq,YRi),YRi),x),x):

if has(deg1,degset) then deg1:=0 fi:

if has(deg2,degset) then deg2:=0 fi:

if has(deg3,degset) then deg3:=0 fi:

deg1-max(deg2,deg3):

end:

#-----------------------------------------------------------------------------

`MultInt/spec`:=proc(a,degs,tdeg,davar)

local i,j,gu,indxlst:

gu:=davar:

if nops(degs)=1 then

 for i from 0 to degs[1] do

  gu:=subs(cat(a,i)=0,gu):

 od:

else

 indxlst:=`MultInt/makindx`(degs,tdeg):

 for i from 1 to nops(indxlst) do

    gu:=subs((cat(a,(cat(seq(indxlst[i][j],j=1..nops(indxlst[i]))))))=0,gu):

 od:

fi:

 RETURN(gu):

end:

#------------------------------------------------------------------------------
# `MultInt/genpol`: generates a generic polynomial with degree degs[i], in intnvars[i]
#         total degree tdeg, the coefficients are a.i.j...

`MultInt/genpol`:=proc(a,vars,degs,tdeg)

local i,j,gu,gu1,indxlst,indx:

gu1:={}:

gu:=0:

if nops(vars)=1 then

 for i from 0 to degs[1] do

    gu:=gu+(cat(a,i))*vars[1]^i:

    gu1:=gu1 union {cat(a,i)}:

 od:

else
 
 indxlst:=`MultInt/makindx`(degs,tdeg):

 for i from 1 to nops(indxlst) do

  indx:=(cat(seq(op(j,indxlst[i]),j=1..nops(indxlst[i])))):

  j:='j':

  gu:=gu + (cat(a,(indx)))*product(vars[j]^op(j,indxlst[i]),j=1..nops(vars)):

  gu1:=gu1 union {cat(a,(indx))}:

 od:

 fi:
 
 RETURN(gu,gu1):

end:

#------------------------------------------------------------------------------
# `MultInt/ptorpols1`: given a differential equation in Ri(y1,...), YRi:=D_yi Ri, 
#          and prescribed degrees degR(list of degrees w.r.t each var yi)
#          and prescribed mainorder, solves a given equation e=0.

`MultInt/ptorpols1`:=proc(e,intnvars,degR,tdegR,mainorder,N)

local i,j,eq,equ,vars1,vars,poly,polylst,Ypolylst,lu,ope:

ope:=1:   

vars1:={}:

for i from 1 to mainorder do

 ope:=ope+(cat(_a_,i))*N^i:

 vars1:=vars1 union {cat(_a_,i)}:

od:

vars:=vars1:

eq:=e:

polylst:=[]:

for i from 1 to nops(intnvars) do

 poly:=`MultInt/genpol`(cat(_A_,i),intnvars,degR[i],tdegR[i]):

 polylst:=[op(polylst),poly[1]]:

 vars:=vars union poly[2]:

od:

Ypolylst:=[]:

for i from 1 to nops(intnvars) do

 Ypolylst :=[op(Ypolylst),diff(polylst[i],intnvars[i])]:

od:

for i from 1 to nops(polylst) do

 eq:=subs((cat(_R_,i)) = polylst[i],eq):

od:

for i from 1 to nops(Ypolylst) do

 eq:=subs((cat(YR_,i)) = Ypolylst[i],eq):

od: 

equ:={coeffs(collect(eq,intnvars,'distributed'),intnvars)}:
 
print(cat(` ... solving `,nops(equ), ` equations for `, nops(vars),

          ` unknowns`)):

lu:=readlib(SolveTools:-Linear)(equ,vars):

if lu=NULL then

  RETURN(0):

fi:

polylst:=subs(lu,polylst):

ope:=subs(lu,ope):

for i from 1 to nops(polylst) do

 for j from 1 to nops(intnvars) do

  polylst:=subsop(i=`MultInt/spec`(cat(_A_,j),degR[j],tdegR[j],polylst[i]),polylst):

 od:

od:


for j from 1 to nops(intnvars) do

   ope:=`MultInt/spec`(cat(_A_,j),degR[j],tdegR[j],ope):
  
od:

if has(ope,vars1) then

 ope:=subs({seq(cat(_a_,j) = 0, j=1..mainorder)}, ope):
 
fi:

for i from 1 to nops(polylst) do

  polylst:=subsop(i=factor(polylst[i]),polylst):

od:

RETURN([ope,polylst]):

end:

#-----------------------------------------------------------------------------
# `MultInt/ptorpolstik1`, like ptorpols but with random assignments of the parameters in
# the list axvars

`MultInt/ptorpolstik1`:=proc(e,intnvars,degR,tdegR,mainorder, axvars)

local i,eq,equ,vars,poly,polylst,Ypolylst,lu,ra:

ra:=2:

eq:=e:

for i from 1 to nops(axvars) do

 eq:=subs(op(i,axvars)=ra,eq):

od:

vars:={}:

for i from 1 to mainorder do

 vars:=vars union {cat(_a_,i)}:

od:

polylst:=[]:

for i from 1 to nops(intnvars) do

 poly:=`MultInt/genpol`(cat(_A_,i),intnvars,degR[i],tdegR[i]):

 polylst:=[op(polylst),poly[1]]:

 vars:=vars union poly[2]:

od:

Ypolylst:=[]:

for i from 1 to nops(intnvars) do

 Ypolylst :=[op(Ypolylst),diff(polylst[i],intnvars[i])]:

od:

for i from 1 to nops(polylst) do

 eq:=subs((cat(_R_,i)) = polylst[i],eq):

od:

for i from 1 to nops(Ypolylst) do

 eq:=subs((cat(YR_,i))=Ypolylst[i],eq):

od:

equ:={coeffs(collect(eq,intnvars,'distributed'),intnvars)}:

lu:=readlib(SolveTools:-Linear)(equ,vars):

if lu<>NULL then

  RETURN(1):

else

  RETURN(0):

fi:

end:

#------------------------------------------------------------------------------
# `MultInt/ptorpols2`: given a differential equation in Ri(y1,...), YRi:=D_yi Ri, 
#          and prescribed degrees degR(list of degrees w.r.t each var yi)
#          and prescribed mainorder, solves a given equation e=0.

`MultInt/ptorpols2`:=proc(e,intnvars,degR,tdegR,mainords_or_strlst)

local i,j,k,eq,equ,indxlst,vars1,vars,poly,polylst,Ypolylst,lu,ope,indx,cfindx:

ope:=1:   

vars1:={}:

if type(mainords_or_strlst,list(nonnegint)) and nops(mainords_or_strlst)>0 then

 indxlst:=`MultInt/makindx2`(mainords_or_strlst):
 
else

 indxlst:=mainords_or_strlst:
 
fi:

for i from 1 to nops(indxlst) do

 indx:= op(i,indxlst):
 
 cfindx:=(cat(seq(op(j,indx),j=1..nops(indx)))): 
 
 ope:=ope + (cat(_a_,cfindx))*product('(cat(_N_,k))^op(k,indx)','k'=1..nops(indx)):
 
 vars1:=vars1 union {cat(_a_,cfindx)}:
 
od:

vars:=vars1:

eq:=e:

polylst:=[]:

for i from 1 to nops(intnvars) do

 poly:=`MultInt/genpol`(cat(_A_,i),intnvars,degR[i],tdegR[i]):

 polylst:=[op(polylst),poly[1]]:

 vars:=vars union poly[2]:

od:

Ypolylst:=[]:

for i from 1 to nops(intnvars) do

 Ypolylst :=[op(Ypolylst),diff(polylst[i],intnvars[i])]:

od:

for i from 1 to nops(polylst) do

 eq:=subs((cat(_R_,i)) = polylst[i],eq):

od:

for i from 1 to nops(Ypolylst) do

 eq:=subs((cat(YR_,i)) = Ypolylst[i],eq):

od: 

equ:={coeffs(collect(eq,intnvars,'distributed'),intnvars)}:
 
print(cat( `  ... solving  `,nops(equ),  ` equations for  `, nops(vars),

           ` unknowns `)):


lu:=readlib(SolveTools:-Linear)(equ,vars):

if lu=NULL then

  RETURN(0):

fi:

polylst:=subs(lu,polylst):


ope:=subs(lu,ope):

for i from 1 to nops(polylst) do

 for j from 1 to nops(intnvars) do

  polylst:=subsop(i=`MultInt/spec`(cat(_A_,j),degR[j],tdegR[j],polylst[i]),polylst):

 od:

od:


for j from 1 to nops(intnvars) do

  ope:=`MultInt/spec`(cat(_A_,j),degR[j],tdegR[j],ope):

od:

if has(ope,vars1) then

 ope:=subs({seq(op(j,vars1) = 0, j=1..nops(vars1))}, ope):

fi:

polylst:=subs({seq(op(j,vars1) = 0, j=1..nops(vars1))}, polylst):

for i from 1 to nops(polylst) do

  polylst:=subsop(i=factor(polylst[i]),polylst):

od:

RETURN([ope,polylst]):

end:


#-----------------------------------------------------------------------------

# `MultInt/ptorpolstik2`, like ptorplos2 but with random assignments of 
#  the parameters in the list axvars

`MultInt/ptorpolstik2`:=proc(e,intnvars,degR,tdegR,mainords_or_strlst, axvars)

local i,j,eq,equ,vars,poly,polylst,Ypolylst,lu,ra,indxlst,indx,cfindx:

ra:=2:

eq:=e:

for i from 1 to nops(axvars) do

 eq:=subs(op(i,axvars)=ra,eq):

od:

vars:={}:


if type(mainords_or_strlst,list(nonnegint)) and nops(mainords_or_strlst)>0 then

 indxlst:=`MultInt/makindx2`(mainords_or_strlst):
 
else

 indxlst:=mainords_or_strlst:
 
fi:

for i from 1 to nops(indxlst) do

 indx:= op(i,indxlst):
 
 cfindx:=(cat(seq(op(j,indx),j=1..nops(indx)))): 
 
 vars:=vars union {cat(_a_,cfindx)}:

od:

polylst:=[]:

for i from 1 to nops(intnvars) do

 poly:=`MultInt/genpol`(cat(_A_,i),intnvars,degR[i],tdegR[i]):

 polylst:=[op(polylst),poly[1]]:

 vars:=vars union poly[2]:

od:

Ypolylst:=[]:

for i from 1 to nops(intnvars) do

 Ypolylst :=[op(Ypolylst),diff(polylst[i],intnvars[i])]:

od:

for i from 1 to nops(polylst) do

 eq:=subs((cat(_R_,i)) = polylst[i],eq):

od:

for i from 1 to nops(Ypolylst) do

 eq:=subs((cat(YR_,i))=Ypolylst[i],eq):

od:

equ:={coeffs(collect(eq,intnvars,'distributed'),intnvars)}:

lu:=readlib(SolveTools:-Linear)(equ,vars):

if lu <> NULL then

  RETURN(1):
   
else

  RETURN(0):

fi:

end:



#------------------------------------------------------------------------------

`MultInt/symR`:=proc(Rlst,intnvars)

local i, j, k, f, A, B,fA,fAj,Sn, Lst,pos, SR,Rpos, nv,indx:

with(combinat):

nv := nops(intnvars):

A:=[seq([j,seq(i,i=1..nv)],j=1..nv)]:

Sn := permute(nv): Lst:=[]:

for i from 1 to nops(Sn) do

 f:=op(i,Sn):
 
 f:={seq(j=f[j],j=1..nops(f))}:
 
 fA:=subs(f,A):
 
 Lst:=[op(Lst),fA]:
  
od:

B:=[]: 

for i from 1 to nops(Lst) do

 fA:=op(i,Lst):
 
 for j from 1 to nops(fA) do
  
  fAj:=fA[j]:
 
  if op(1,fAj)=1 then 
  
   indx:=subsop(1=j,fAj):
  
   B:=[op(B),indx]:
   
  fi:
 
 od:

od:
 
SR:=0:

for i from 1 to nops(B) do
 
  indx := op(i,B): 
  
  pos:=op(1,indx): indx:=subsop(1=NULL,indx):
   
  Rpos:=subs({seq(op(k,intnvars)=op(op(k,indx),intnvars),k=1..nv)},
           
           op(pos,Rlst)):
           
  SR:=SR + Rpos:
  
od:

factor(SR/nv!):

end:


#-------------------------------------------------------------------------------
`MultInt/findrec10`:=proc(p,f,mainvar,intnvars,axvars,mainorder,polylst,fnam,st)

local eq,gu,su,i,j,degR,tdegR,Rlst,SR,ope,mekh,N,prts,count, fvars, rhside, deltintnvar:

eq:=`MultInt/fa1`(p,f,mainvar,intnvars,mainorder):

su:=[seq(0,i=1..nops(intnvars))]:

for i from 1 to nops(polylst) do

  su:=subsop(i=1/polylst[i],su):

od:

eq:=`MultInt/cv`(su,eq,intnvars):

if type(eq,`+`) then

 prts:=[op(eq)]:

 eq:=0:

 for i from 1 to nops(prts) do

  eq:=factor(eq +  factor(op(i,prts))):
 
 od:
 
fi:

eq:=numer(eq):

degR:=[seq([seq(0,i=1..nops(intnvars))],j=1..nops(intnvars))]:

tdegR:=[seq(0,i=1..nops(intnvars))]:

for i from 1 to nops(intnvars) do

 
degR:=subsop(i=[seq(`MultInt/finddeg`(eq,cat(_R_,i),cat(YR_,i),op(j,intnvars)),

              j=1..nops(intnvars))],degR):
 
 tdegR:=subsop(i=convert(op(i,degR),`+`),tdegR):
 
od:

#print(``): 

#lprint(` ... checking existence of solution by random substitution for auxiliary vars`):

count:=1: 

if mainorder=0 then count:=4: fi:

gu:=`MultInt/ptorpolstik1`(eq,intnvars,degR,tdegR, mainorder, axvars):

while  gu=0 and count <= 4 do

 for i from 1 to nops(intnvars) do

  degR:=subsop(i=[seq(op(j,op(i,degR)) + 1,j=1..nops(intnvars))],degR):

  tdegR:=subsop(i=op(i,tdegR) + nops(intnvars),tdegR):

 od:

 gu:=`MultInt/ptorpolstik1`(eq,intnvars,degR,tdegR,mainorder,axvars):
 
 count:=count+1:
 
od:

if gu=0 then

  RETURN(NULL):

fi:


 
#Now trying to get a non-zero annihilator ...

gu:=`MultInt/ptorpols1`(eq,intnvars,degR,tdegR,mainorder,N);

if gu=0 then

 for i from 1 to nops(intnvars) do

  degR:=subsop(i=[seq(op(j,op(i,degR))+1,j=1..nops(intnvars))],degR):

  tdegR:=subsop(i=op(i,tdegR)+nops(intnvars),tdegR):

 od:

 gu:=`MultInt/ptorpolstik1`(eq,intnvars,degR,tdegR, mainorder,axvars):

 if gu=0 then

  RETURN(NULL):

 else

  gu:=`MultInt/ptorpols1`(eq,intnvars,degR,tdegR,mainorder,N):

  if gu=0 then
    
   RETURN(NULL):

  fi:

 fi:

fi:


if gu=0 then

 RETURN(NULL):

fi:

if gu = 1 then

  gu:=`MultInt/ptorpols1`(eq,intnvars,degR,tdegR,mainorder,N):

  if gu=0 then
    
    RETURN(NULL):

   fi:
fi:

ope:=op(1,gu):

Rlst:=[]:

for i from 1 to nops(intnvars) do

 Rlst:=[op(Rlst), op(i,op(2,gu))/op(i,polylst)]:

od:

ope:=normal(ope):

mekh:=denom(ope):

ope:=numer(ope):

ope:=`MultInt/pashet1`(ope,N,mainvar,intnvars,fnam):

for i from 1 to nops(intnvars) do

 Rlst:=subsop(i=factor(op(i,Rlst)*mekh/p), Rlst):

od:

if nops(intnvars) >1 and  `MultInt/sym`(p*f,intnvars) then

 SR := `MultInt/symR`(Rlst,intnvars):
 
  printf(`%s%s%.2f%s\n`,` Cheers! for the success. `,  ` CPU Time : `,
            
             time() - st, ` seconds.`): print(``):  
           

fvars:=op(intnvars):
rhside:=0: 
for i from 1 to nops(intnvars) do 
deltintnvar:=op(subsop(i=NULL,intnvars)):
rhside:=rhside + (Diff(_R(intnvars[i],deltintnvar)*fnam(mainvar,fvars),intnvars[i])):
od:             
             
RETURN([ope=rhside,_R(fvars) = SR]):



#RETURN([ope=sum('Diff(_R(intnvars[i],op(subsop(i=NULL,intnvars)))*

#           fnam(mainvar,op(intnvars)),intnvars[i])','i'=1..nops(intnvars)),

#           _R(op(intnvars)) = SR]):

fi:
 
printf(`%s%s%.2f%s\n`,` Cheers! for the success. `,  ` CPU Time : `,
            
             time() -st, ` seconds.`): print(``):

fvars:=op(intnvars):
rhside:=0: 
for i from 1 to nops(intnvars) do 
rhside:=rhside + (Diff(Rlst[i]*fnam(mainvar, fvars),intnvars[i])):
od:             
             
RETURN(ope=rhside):

end:


#-----------------------------------------------------------------------------

`MultInt/findrec101`:=proc(p,f,mainvar,intnvars,axvars,mainorder,fnam,st)

local i,j,den,den1,den11,den12,den2,dfbyf,denlst,result,polylst,maxd,factr,factrs,

      factrs2,indxlst,indx:

result:=NULL: den:=`MultInt/extrctdenom`(f,intnvars):

factrs:=op(2,factors(den)):

maxd:=0:

for i from 1 to nops(factrs) do
 
 factr:=collect(op(1,op(i,factrs))^op(2,op(i,factrs)),intnvars,'distributed'):

 maxd:=max(maxd,degree(factr,{op(intnvars)})):
 
od:


if mainorder=0  then maxd:=0 fi:

 if not type(den,`+`) then 

  den1:=select(type,den,{name ^ nonnegint, name ^ 1}):

 else

  den1:=den:
 
 fi:

 if  result=NULL and den1<>NULL and has(den1,intnvars) then

for i from 0 to mainorder + maxd while result=NULL do

 den11:=factor(den1^i):
 
 polylst:=[seq(den11,j=1..nops(intnvars))]:


result:=`MultInt/findrec10`(p,f,mainvar,intnvars,axvars,mainorder,polylst,fnam,st):

od:

fi:

if result=NULL and den1<>NULL and has(den1,intnvars) then

  denlst:=[]:

 for i from 1 to nops(intnvars) do

   dfbyf:=diff(simplify(log(f)),intnvars[i]):  

   dfbyf:=normal(dfbyf):

   denlst:=[op(denlst),denom(dfbyf)]:

 od:

 polylst:=[seq(normal(op(i,denlst)/gcd(den1,op(i,denlst))),i=1..nops(denlst))]:
   

result:=`MultInt/findrec10`(p,f,mainvar,intnvars,axvars,mainorder,polylst,fnam,st):

fi:
 
 
if result=NULL and den1<>NULL and has(den1,intnvars) then

 den2:=factor(den/den1):

 factrs:=op(2,factors(den2)): factrs2:=[]:

 for i from 1 to nops(factrs) while result=NULL do

  factr:=op(i,factrs):
 
  factr:=op(1,factr)^(op(2,factr)): factr:=factor(factr):

  if has(factr,intnvars) then

   factrs2:=[op(factrs2),factr]:

   polylst:=[seq(factr,j=1..nops(intnvars))]:
 
  
result:=`MultInt/findrec10`(p,f,mainvar,intnvars,axvars,mainorder,polylst,fnam,st):  
  fi:
 
 od:

 if nops(factrs2)>1 and result=NULL then
  
  with(combinat):
   
  indxlst:=choose(nops(factrs2),2):
  
  for i from 1 to nops(indxlst) while result=NULL do
   
   indx:=op(i,indxlst):
    
   factr:= op(indx[1],factrs2)*op(indx[2],factrs2):
   
   factr:=factor(factr):
    
   polylst:=[seq(factr,j=1..nops(intnvars))]:
  
result:=`MultInt/findrec10`(p,f,mainvar,intnvars,axvars,mainorder,polylst,fnam,st): 
   od:
   
 fi:

fi:
 
if result=NULL and den1<>NULL and has(den1,intnvars) then

  for i from 0 to mainorder + maxd while result=NULL do

    den12:=factor(den1^i*den2):
   
    polylst:=[seq(den12,j=1..nops(intnvars))]:

   
result:=`MultInt/findrec10`(p,f,mainvar,intnvars,axvars,mainorder,polylst,fnam,st):

  od:

fi:

if result=NULL and den1<>NULL and has(den1,intnvars) then

  den2:=factor(den/den1):
 
  for i from 2 to mainorder + maxd while result=NULL do
  
   den12:=factor(den1*den2^i):
   
   polylst:=[seq(den12,j=1..nops(intnvars))]:

  result:=`MultInt/findrec10`(p,f,mainvar,intnvars,axvars,mainorder,polylst,fnam,st):

  od: 

fi:

if result=NULL then
 
  for i from 0 to mainorder+maxd  while result=NULL do

   den1:=factor(den^i):
   
   polylst:=[seq(den1, j=1..nops(intnvars))]:
  
result:=`MultInt/findrec10`(p,f,mainvar,intnvars,axvars,mainorder,polylst,fnam,st):

  od:
   
fi:

if result=NULL then
  
 print(cat(`could not find a non-zero recurrence-differential(WZ)`,

` equation with order `, mainorder )): 

fi:

result:
 
end:  


#-----------------------------------------------------------------------------

`MultInt/findrec102`:=proc(p,f,mainvar,intnvars,axvars,polylst,st)

local result,maxorder,mainorder:

maxorder:=6:  result:=NULL:

for mainorder from 0 to maxorder while result=NULL do

 print(cat(` ... trying to find a non-zero recurrence-differential(WZ)`,

` equation with order `, mainorder)): 
 
result:=`MultInt/findrec10`(p,f,mainvar,intnvars,axvars,mainorder,polylst,fnam,st):

od:

if result=NULL then
  
  print(cat(` could not get a non-zero recurrence eq with max order `,
  
         maxorder)): 
fi:

result:

end:  

#-----------------------------------------------------------------------------

`MultInt/findrec103`:=proc(p,f,mainvar,intnvars,axvars,fnam,st)

local mainorder,maxorder,result:

maxorder:=6: result:=NULL:

for mainorder from 0 to maxorder while result=NULL do

 print(cat(`... trying to find a non-zero recurrence-differential(WZ)`,

` equation with order `, mainorder)): 
          
 result:=`MultInt/findrec101`(p,f,mainvar,intnvars,axvars,mainorder,fnam,st):

od:

if result=NULL then
  
  print(cat(` could not get a non-zero recurrence-differential(WZ)`,

 ` equation with max. order `,  maxorder )): 

fi:

result:
 
end:  

#-------------------------------------------------------------------------------
`MultInt/findrec20`:=proc(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st)

local eq,gu,su,i,j,degR,tdegR,Rlst,SR,ope,mekh,prts,count,mvars, fvars, rhside, deltintnvar:

eq:=`MultInt/fa2`(p,f,mainvars,intnvars,mainorders):

su:=[seq(0,i=1..nops(intnvars))]:

for i from 1 to nops(polylst) do

  su:=subsop(i=1/polylst[i],su):

od:

eq:=`MultInt/cv`(su,eq,intnvars):

if type(eq,`+`) then

 prts:=[op(eq)]:

 eq:=0:

 for i from 1 to nops(prts) do

  eq:=factor(eq +  factor(op(i,prts))):

 od:
 
fi:

eq:=numer(eq):

degR:=[seq([seq(0,i=1..nops(intnvars))],j=1..nops(intnvars))]:

tdegR:=[seq(0,i=1..nops(intnvars))]:

for i from 1 to nops(intnvars) do

 
degR:=subsop(i=[seq(`MultInt/finddeg`(eq,cat(_R_,i),cat(YR_,i),op(j,intnvars)),

             j=1..nops(intnvars))],degR):
 
 tdegR:=subsop(i=convert(op(i,degR),`+`),tdegR):
 
od:

#print(``): 

#lprint(` ... checking existence of solution by random substitution for 
#        the auxiliary vars`):

count:=1: 

if nops({op(mainorders)})=1 and op({op(mainorders)})=0 then

 count:=4:
 

fi:

gu:=`MultInt/ptorpolstik2`(eq,intnvars,degR,tdegR, mainorders, axvars):

while  gu=0 and count<=4 do

 for i from 1 to nops(intnvars) do

  degR:=subsop(i=[seq(op(j,op(i,degR)) + 1,j=1..nops(intnvars))],degR):

  tdegR:=subsop(i=op(i,tdegR) + nops(intnvars),tdegR):

 od:

 gu:=`MultInt/ptorpolstik2`(eq,intnvars,degR,tdegR,mainorders,axvars):
 
 count:=count + 1:
 
od:

if gu=0 then

 RETURN(NULL):

fi:


#High probability! (but no guarantee) to get a non-zero annihilator. 
 
#Now trying to get a non-zero annihilator ...

 gu:=`MultInt/ptorpols2`(eq,intnvars,degR,tdegR,mainorders);

 if gu=0 then

  for i from 1 to nops(intnvars) do

   degR:=subsop(i=[seq(op(j,op(i,degR))+1,j=1..nops(intnvars))],degR):

   tdegR:=subsop(i=op(i,tdegR)+nops(intnvars),tdegR):

  od:

   gu:=`MultInt/ptorpolstik2`(eq,intnvars,degR,tdegR, mainorders,axvars):

   if gu=0 then

    RETURN(NULL):

  else

   gu:=`MultInt/ptorpols2`(eq,intnvars,degR,tdegR,mainorders):

   if gu=0 then
    
    RETURN(NULL):

   fi:

  fi:

 fi:

if gu=0 then

 RETURN(NULL):

fi:

if gu=1 then

 gu:=`MultInt/ptorpols2`(eq,intnvars,degR,tdegR,mainorders):

 if gu=0 then RETURN(NULL) fi:

fi:

ope:=gu[1]:

Rlst:=[]:

for i from 1 to nops(intnvars) do

 Rlst:=[op(Rlst), op(i,op(2,gu))/op(i,polylst)]:

od:

ope:=normal(ope):

mekh:=denom(ope):

ope:=numer(ope):

ope:=`MultInt/pashet2`(ope,mainvars,intnvars,fnam):

for i from 1 to nops(intnvars) do

 Rlst:=subsop(i=factor(Rlst[i]*mekh/p), Rlst):

od:

if nops(intnvars) > 1 and  `MultInt/sym`(p*f,intnvars) then

 SR := `MultInt/symR`(Rlst,intnvars):
 
 printf(`%s%s%.2f%s\n`,` Cheers! for the success. `,  ` CPU Time : `,
            
             time() -st, ` seconds.`):  print(``):
             
             
             
fvars:=op(intnvars):
rhside:=0: 
for i from 1 to nops(intnvars) do 
deltintnvar:=op(subsop(i=NULL,intnvars)):
rhside:=rhside + (Diff(_R(intnvars[i],deltintnvar)*fnam(mainvar,fvars),intnvars[i])):
od:             
             
RETURN([ope=rhside,_R(fvars) = SR]):

 
#RETURN([ope=sum('Diff(_R(intnvars[i],op(subsop(i=NULL,intnvars)))*

 #            fnam(op(mainvars),op(intnvars)),intnvars[i])','i'=1..nops(intnvars)),

  #           _R(op(intnvars)) = SR]):

fi:
 
printf(`%s%s%.2f%s\n`,` Cheers! for the success. `,  ` CPU Time: `,
            
             time() - st, ` seconds.`):  

fvars:=op(intnvars):
mvars:=op(mainvars):
rhside:=0: 
for i from 1 to nops(intnvars) do 
deltintnvar:=op(subsop(i=NULL,intnvars)):
rhside:=rhside + Diff(Rlst[i]*fnam(mvars, fvars),intnvars[i]):
od:             
   
RETURN(ope=rhside):

#RETURN(ope=sum('Diff(Rlst[i]*fnam(op(mainvars),op(intnvars)),intnvars[i])',
#      'i'=1..nops(intnvars))):

end:


#-----------------------------------------------------------------------------

`MultInt/findrec201`:=proc(p,f,mainvars,intnvars,axvars,mainorders,fnam,st)

local i,j,den,den1,den11,den12,den2,result,polylst,maxd,maxord,factr,factrs,factrs2,

      indxlst, dfbyf, indx,denlst:

result:=NULL: den:=`MultInt/extrctdenom`(f,intnvars):

factrs:=op(2,factors(den)):

maxd:=0:

for i from 1 to nops(factrs) do

 factr:=collect(op(1,op(i,factrs))^op(2,op(i,factrs)),intnvars,'distributed'):

 maxd:=max(maxd,degree(factr,{op(intnvars)})):
 
od:

maxord:=max(op(mainorders)):

if maxord=0  then maxd:=0 fi:

den1:=select(type,den,{name ^ nonnegint, name ^ 1}):

if  result=NULL and den1<>NULL and has(den1,intnvars) then

 for i from 0 to maxord + maxd while result=NULL do

  den11:=factor(den1^i):
 
  polylst:=[seq(den11,j=1..nops(intnvars))]:

 result:=`MultInt/findrec20`(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st):

 od:

fi:

if result=NULL and den1<>NULL and has(den1,intnvars) then

 denlst:=[]:

 for i from 1 to nops(intnvars) do

  dfbyf:=diff(simplify(log(f)),intnvars[i]):  

  dfbyf:=normal(dfbyf):

  denlst:=[op(denlst),denom(dfbyf)]:

 od:

 polylst:=[seq(normal(op(i,denlst)/gcd(den1,op(i,denlst))),i=1..nops(denlst))]:
   
result:=`MultInt/findrec20`(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st):

fi: 

if result=NULL and den1<>NULL and has(den1,intnvars) then

 den2:=factor(den/den1):

 factrs:=op(2,factors(den2)): factrs2:=[]:

 for i from 1 to nops(factrs) while result=NULL do

  factr:=op(i,factrs):
 
  factr:=op(1,factr)^(op(2,factr)): factr:=factor(factr):

  if has(factr,intnvars) then

   factrs2:=[op(factrs2),factr]:
 
   polylst:=[seq(factr,j=1..nops(intnvars))]:
 
  result:=`MultInt/findrec20`(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st):

  fi:
 
 od:

 if nops(factrs2)>1 and result=NULL then
  
   with(combinat):
   
  indxlst:=choose(nops(factrs2),2):
  
  for i from 1 to nops(indxlst) while result=NULL do
   
   indx:=op(i,indxlst):
    
   factr:= op(indx[1],factrs2)*op(indx[2],factrs2):
   
   factr:=factor(factr):
    
   polylst:=[seq(factr,j=1..nops(intnvars))]:

  result:=`MultInt/findrec20`(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st):

  od:
   
 fi:
 
 for i from 0 to maxord + maxd while result=NULL do
     
   den12:=factor(den1^i*den2):
   
   polylst:=[seq(den12,j=1..nops(intnvars))]:

result:=`MultInt/findrec20`(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st):

  od:

fi:

if result=NULL and den1<>NULL and has(den1,intnvars) then

  den2:=factor(den/den1):

  for i from 2 to maxord + maxd while result=NULL do
  
   den12:=factor(den1*den2^i):
   
   polylst:=[seq(den12,j=1..nops(intnvars))]:

  result:=`MultInt/findrec20`(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st):

  od: 

fi:

if result=NULL then
 
  for i from 0 to maxord + maxd while result=NULL do

   den1:=factor(den^i):
   
   polylst:=[seq(den1,j=1..nops(intnvars))]:

  result:=`MultInt/findrec20`(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st):

  od:
   
fi:

if result=NULL then

 if nops(mainvars)=1 then  

  print(cat(`could not find a non-zero recurrence-differential(WZ)`,
 ` equation with order `,convert(op(mainorders),string), ` in `,convert(op(mainvars),string))):
         
 else
 
  print(cat(`could not find a non-zero recurrence-differential(WZ)`,
` equation with order `, convert(mainorders,string), ` in `,convert(mainvars,string))):      
 fi:
          
fi:

result:

end:  


#-----------------------------------------------------------------------------

`MultInt/orderindx`:=proc(indx,nmainv,n)

local indx1,newindx,i,j,temp,indx2:

indx1:=indx: newindx:=[]:

for i from 1 to nmainv do

 temp:=[]:
  
 for j from 1 to nops(indx1)  do
   
  if op(i,indx1[j]) = n - 1 then
  
   indx2:=subsop(i=n,indx1[j]):
   
   newindx:= [op(newindx),indx2]:
   
   temp:=[op(temp),indx2]:
      
  fi:
  
 od:
 
 indx1:=temp:
 
od:

newindx:

end:
 
#-------------------------------------------------------------------------------

`MultInt/findrec202`:=proc(p,f,mainvars,intnvars,axvars,polylst,fnam,st)

local result,maxorder,mainorders,indx,i:

maxorder:=6:  indx:=[]:

mainorders:=[seq(0,i=1..nops(mainvars))]: 

indx:=[op(indx),mainorders]:

for i from 1 to nops(mainvars) do

  mainorders:=subsop(i=1,mainorders):
 
 indx:=[op(indx),mainorders]:
 
od:

result:=`MultInt/findrec2021`(p,f,mainvars,intnvars,axvars,polylst,indx,fnam,st):

if result=NULL then

  for i from 2 to maxorder while result=NULL do
  
   indx:=`MultInt/orderindx`(indx,nops(mainvars),i):

result:=`MultInt/findrec2021`(p,f,mainvars,intnvars,axvars,polylst,indx,fnam,st):   

  od:

fi:

result:

end:  

#-----------------------------------------------------------------------------
`MultInt/findrec2021`:=proc(p,f,mainvars,intnvars,axvars,polylst,indx,fnam,st)

local result,mainorders,morders,i,j:

result:=NULL:

if not `MultInt/sym`(p*f,mainvars) then

 with(combinat):
 
 for i from 1 to nops(indx)  while result=NULL do

   morders:=permute(op(i,indx)):
  
   for j from 1 to nops(morders) while result=NULL do
  
    mainorders:=op(j,morders):
    
    if nops(mainvars)=1 then
  
      print(cat(`... trying to get anon-zero recurrence-differential(WZ)`,

      ` equation with order `, convert(op(mainorders),string), ` in `,      

              convert(op(mainvars),string))):
         
    else
    
      print(cat(` ... trying to get non-zero a non-zero recurrence-differential(WZ)`,
  
     ` equation with order `, convert(mainorders,string), ` in `, 
         
         convert(mainvars,string))):
         
    fi:
       
result:=`MultInt/findrec20`(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st):
    
  od:
  
 od:
  
else
 
  for i from 1 to nops(indx)  while result=NULL do

    mainorders:=op(i,indx):
    
    if nops(mainvars)=1 then
  
     print(cat(` ... trying to get a non-zero recurrence-differential(WZ)`,

     ` equation with order `,convert(op(mainorders),string), ` in `, 

       convert(op(mainvars),string))):
         
    else
    
      print(cat(` ... trying to get a non-zero recurrence-differential(WZ)`,

      ` equation with order `,convert(mainorders,string), ` in `, 
         
       convert(mainvars,string))):
         
    fi:
    
    
result:=`MultInt/findrec20`(p,f,mainvars,intnvars,axvars,mainorders,polylst,fnam,st):
   
  od:
  
fi:

result:

end:

#-----------------------------------------------------------------------------

`MultInt/findrec203`:=proc(p,f,mainvars,intnvars,axvars,fnam,st)

local mainorders,maxorder,result,indx,i:

maxorder:=6: indx:=[]:

mainorders:=[seq(0,i=1..nops(mainvars))]: 

indx:=[op(indx),mainorders]:

for i from 1 to nops(mainvars) do

 mainorders:=subsop(i=1,mainorders):
 
 indx:=[op(indx),mainorders]:
 
od:

result:=`MultInt/findrec2031`(p,f,mainvars,intnvars,axvars,indx,fnam,st):

if result=NULL then

  for i from 2 to maxorder while result=NULL do
  
   indx:=`MultInt/orderindx`(indx,nops(mainvars),i):
   
   result:=`MultInt/findrec2031`(p,f,mainvars,intnvars,axvars,indx,fnam,st):
   
  od:

fi:

result:

end:

#-----------------------------------------------------------------------------
`MultInt/findrec2031`:=proc(p,f,mainvars,intnvars,axvars,indx,fnam,st)

local result,mainorders,morders,i,j:

result:=NULL:

if not `MultInt/sym`(p*f,mainvars) then

 with(combinat):
 
 for i from 1 to nops(indx)  while result=NULL do

   morders:=permute(op(i,indx)):
  
   for j from 1 to nops(morders) while result=NULL do
  
    mainorders:=op(j,morders):
    
    if nops(mainvars)=1 then
  
     print(cat(` ... trying to get a recurrence eq with order `,
  
         convert(op(mainorders),string), ` in `, 
         
         convert(op(mainvars),string))):
         
    else
    
      print(cat(` ... trying to get a recurrence eq with order `,
  
         convert(mainorders,string), ` in `, 
         
         convert(mainvars,string))):
         
    fi:
    
   result:=`MultInt/findrec201`(p,f,mainvars,intnvars,axvars,mainorders,fnam,st): 

    od:
  
 od:
  
else
 
  for i from 1 to nops(indx)  while result=NULL do

    mainorders:=op(i,indx):
    
    if nops(mainvars)=1 then
  
     print(cat(` ... trying to get a non-zero recurrence-differential(WZ)`,

` equation with order `, convert(op(mainorders),string), ` in `, 
         
         convert(op(mainvars),string))):
         
    else
    
      print(cat(` ... trying to get a non-zero recurrence-differential(WZ) `,
  
        ` equation with order `, convert(mainorders,string), ` in `, 
         
         convert(mainvars,string))):
         
    fi:
    
   result:=`MultInt/findrec201`(p,f,mainvars,intnvars,axvars,mainorders,fnam,st): 

   od:
  
fi:

result:

end:

#-------------------------------------------------------------------------------
#`MultInt/hypertorecdiff`
#MultInt[hypertorecdiff]

`MultInt/hypertorecdiff`:=proc()

local
args1,axvars,vars,res,i,p,f,fvar,_fnam,t11,t12,t13,t14,t21,t22,t23,t24,t31,t32,st:

st:=time():

if nargs < 3 then ERROR(` Wrong number of arguments`) fi:

#defining the correct args types

t11:=[algebraic,name,{list(name),function},nonnegint,list(polynom)]:

t12:=[algebraic,name,{list(name),function},nonnegint]:

t13:=[algebraic,name,{list(name),function},list(polynom)]:

t14:=[algebraic,name,{list(name),function}]:


t21:=[algebraic,list(name),{list(name),function}, list(nonnegint),list(polynom)]:

t22:=[algebraic,list(name),{list(name),function}, list(nonnegint)]:

t23:=[algebraic,list(name),{list(name),function}, list(polynom)]:

t24:=[algebraic,list(name),{list(name),function}]:


t31:=[algebraic,list(name),{list(name),function},list(name), list({1,name,`*`,`^`}), list(polynom)]:

t32:=[algebraic,list(name),{list(name),function},list(name),list({1,name,`*`,`^`})]:  if not type([args],{t11,t12,t13,t14,t21,t22,t23,t24,t31,t32}) then  
 ERROR(` Wrong argument types or type mismatch `,args);  fi:  args1:=[args]: fvar:=args1[3]:   if type(fvar,list(name)) then    vars:=fvar: _fnam:=_F:   elif type(fvar,function) and type(op(0,fvar),name) and type([op(fvar)],list(name)) then 
  
 vars:=[op(fvar)]: _fnam:=op(0,fvar):   else      ERROR(cat(` the third argument must be list of variable names  or\n`,         ` a function of variable names `), fvar):  fi:
 if `MultInt/notdistvars`(vars) then     ERROR(cat(` the variable names must be distinct `,vars)):  fi:  if `MultInt/freeintgrand`(op(1,args1),vars) then    ERROR(`The simplified form of  `, op(1,args1),    cat(` does not contain all variables in `,convert(vars,string),`.`)):  fi:  axvars:=`MultInt/extrctauxvars`(op(1,args1),convert(vars,set)):  if type(args1, {t11,t12,t13,t14}) then 

 args1:=subsop(3=op([vars,axvars]),args1):   args1:=[op(args1),_fnam,st]: 
 res:=`MultInt/findrec1`(op(args1)):

elif type(args1, {t21,t22,t23,t24}) then   args1:=subsop(3=op([vars,axvars]),args1): 
 args1:=[op(args1),_fnam,st]:

 res:=`MultInt/findrec2`(op(args1)):

elif type([args],{t31,t32}) then   args1:=subsop(3=op(4,args1),4=op([vars,axvars]),args1): 
  args1:=[op(args1),_fnam,st]:

 res:=`MultInt/findrec3`(op(args1)):

fi:

if res=NULL then RETURN(0) fi:

res:

 
end: 


MultInt[hypertorecdiff]:=proc()

`MultInt/hypertorecdiff`(args):

end:


#-------------------------------------------------------------------------------
#`Mint/findrec1`

`MultInt/findrec1`:=proc()

local args1,nargs1,args2,i,p,f,polylst,t1,t2,t3,t4,result:

#defining the correct args types

t1:=[algebraic,name,list(name),list(name), nonnegint,list(polynom),name,numeric]:

t2:=[algebraic,name,list(name),list(name), nonnegint,name,numeric]:

t3:=[algebraic,name,list(name),list(name), list(polynom),name,numeric]:

t4:=[algebraic,name,list(name),list(name),name,numeric]:

args1:=[args]: nargs1:=nops(args1): 

args2:=subsop(1=NULL,nargs1=NULL,nargs1-1=NULL,args1): 

if type(args1,t1) then   polylst:=`MultInt/checkpoly`(args1[6],nops(args1[3])):   p:=`MultInt/extrctP`(seq(args1[i],i=1..4)):   f:=p[2]: p:=p[1]:

 result:=`MultInt/findrec10`(p,f,op(args2),args1[5],polylst,args1[7],args1[8]):

 if result=NULL then 
 
  print(` could not find the required non-zero rec-diff equation. `):
  
 else 
 
  result:
 
 fi:
 

elif type(args1,t2) then

 
 p:=`MultInt/extrctP`(seq(args1[i],i=1..4)):

 f:=p[2]: p:=p[1]:

 `MultInt/findrec101`(p,f,op(args2),args1[6],args1[7]):

elif type(args1,t3) then

 polylst:=`MultInt/checkpoly`(args2[nops(args2)],nops(args1[3])):

 p:=`MultInt/extrctP`(seq(args1[i],i=1..4)):

 f:=p[2]: p:=p[1]:
 
 `MultInt/findrec102`(p,f,op(args2),polylst,args1[6],args1[7]):

elif type(args1,t4) then

 args2:=subsop(1=NULL,args1):

 p:=`MultInt/extrctP`(seq(args1[i],i=1..4)):

 f:=p[2]: p:=p[1]:
 
 `MultInt/findrec103`(p,f,op(args2),args1[5],args1[6]):

fi:

end:


#--------------------------------------------------------------------------
`MultInt/findrec2`:=proc()

local args1,nargs1,args2,i,p,f,setdf,polylst,t1,t2,t3,t4:

if nargs< 3 then ERROR(` Wrong number of arguments`) fi:

#defining the correct args types
  
t1:=[algebraic,list(name),list(name),list(name),list(nonnegint),list(polynom),name,numeric]:

t2:=[algebraic,list(name),list(name),list(name),list(nonnegint),name,numeric]:

t3:=[algebraic,list(name),list(name),list(name),list(polynom),name,numeric]:

t4:=[algebraic,list(name),list(name),list(name),name,numeric]:

args1:=[args]: nargs1:=nops(args1):

if nops(op(2,args1))=1 and type(args1,t4) then 

 args1:=subsop(2=op(op(2,args1)),args1):
 
 RETURN(`MultInt/findrec1`(op(args1))):

fi:


args2:=subsop(1=NULL,nargs1=NULL,nargs1-1=NULL,args1):

if type(args1,t1) then

  if not `MultInt/equalsize`(args1[2],args1[5]) then 

  ERROR(args1[2], args1[5],` must be of equal size`):

 fi:

 polylst:=`MultInt/checkpoly`(args2[nops(args2)],nops(args1[3])):

 p:=`MultInt/extrctP`(seq(args1[i],i=1..4)):

 f:=p[2]: p:=p[1]:

 `MultInt/findrec20`(p,f,op(args2),args1[5],polylst,args1[7],args1[8]):

elif type(args1,t2) then

 if not `MultInt/equalsize`(args1[2],args1[5]) then 

  ERROR(args1[2], args1[5],` must be of equal size`):

 fi:

 args2:=subsop(1=NULL,args1):

 p:=`MultInt/extrctP`(seq(args1[i],i=1..4)):

 f:=p[2]: p:=p[1]:

 `MultInt/findrec201`(p,f,op(args2),op(nargs1-1,args1),op(nargs1,args1)):

elif type(args1,t3) then

 polylst:=`MultInt/checkpoly`(args2[nops(args2)],nops(args1[3])):

 p:=`MultInt/extrctP`(seq(args1[i],i=1..4)):

 f:=p[2]: p:=p[1]:

 `MultInt/findrec202`(p,f,op(args2),polylst,args1[6],args1[7]):

elif type(args1,t4) then

 args2:=subsop(1=NULL,args1):

 p:=`MultInt/extrctP`(seq(args1[i],i=1..4)):

 f:=p[2]: p:=p[1]:

 `MultInt/findrec203`(p,f,op(args2),args1(nargs1-1),args1[nargs1]):

fi:

end:


#-------------------------------------------------------------------------------
`MultInt/findrec30`:=proc(p,f,mainvars,intnvars,axvars,strlst,polylst,fnam,st)

local eq,gu,su,i,j,degR,tdegR,Rlst,SR,ope,mekh,prts,count,fvars, 
mvars, rhside, deltintnvar:

eq:=`MultInt/fa2`(p,f,mainvars,intnvars,strlst):

su:=[seq(0,i=1..nops(intnvars))]:

for i from 1 to nops(polylst) do

  su:=subsop(i=1/polylst[i],su):

od:

eq:=`MultInt/cv`(su,eq,intnvars):

if type(eq,`+`) then

 prts:=[op(eq)]:

 eq:=0:

 for i from 1 to nops(prts) do

  eq:=factor(eq +  factor(op(i,prts))):

 od:
 
fi:

eq:=numer(eq):

degR:=[seq([seq(0,i=1..nops(intnvars))],j=1..nops(intnvars))]:

tdegR:=[seq(0,i=1..nops(intnvars))]:

for i from 1 to nops(intnvars) do

 degR:=subsop(i=[seq(`MultInt/finddeg`(eq,cat(_R_,i),cat(YR_,i),op(j,intnvars)),
 
               j=1..nops(intnvars))],degR):
 
 tdegR:=subsop(i=convert(op(i,degR),`+`),tdegR):
 
od:


#...checking existence of solution by random substitution for auxiliary vars

count:=1: 

if max(op(map(degree,strlst)))=0 then

 count:=4:
 
fi:

gu:=`MultInt/ptorpolstik2`(eq,intnvars,degR,tdegR, strlst, axvars):

while  gu=0 and count<=4 do

 for i from 1 to nops(intnvars) do

  degR:=subsop(i=[seq(op(j,op(i,degR)) + 1,j=1..nops(intnvars))],degR):

  tdegR:=subsop(i=op(i,tdegR) + nops(intnvars),tdegR):

 od:

 gu:=`MultInt/ptorpolstik2`(eq,intnvars,degR,tdegR,strlst,axvars):
 
 count:=count + 1:
 
od:

if gu=0 then

   RETURN(NULL):

fi:

gu:=`MultInt/ptorpolstik2`(eq,intnvars,degR,tdegR,strlst,axvars):

if gu=0 then

 for i from 1 to nops(intnvars) do

  degR:=subsop(i=[seq(op(j,op(i,degR)) + 1,j=1..nops(intnvars))],degR):

  tdegR:=subsop(i=op(i,tdegR) + nops(intnvars),tdegR):

 od:

 gu:=`MultInt/ptorpolstik2`(eq,intnvars,degR,tdegR,strlst,axvars):

 if gu=0 then

   RETURN(NULL):

 fi:

# High probability! (but no guarantee) to get a non-zero annihilator. 

#Now trying to get a non-zero annihilator ...

 gu:=`MultInt/ptorpols2`(eq,intnvars,degR,tdegR,strlst);

 if gu=0 then

  for i from 1 to nops(intnvars) do

   degR:=subsop(i=[seq(op(j,op(i,degR))+1,j=1..nops(intnvars))],degR):

   tdegR:=subsop(i=op(i,tdegR)+nops(intnvars),tdegR):

  od:

   gu:=`MultInt/ptorpolstik2`(eq,intnvars,degR,tdegR, strlst,axvars):

   if gu=0 then

    RETURN(NULL):

  else

   gu:=`MultInt/ptorpols2`(eq,intnvars,degR,tdegR,strlst):

   if(gu=0) then
    
    RETURN(NULL):

   fi:

  fi:

 fi:

fi:

if gu=0 then

 RETURN(NULL):

fi:

if gu=1 then

 gu:=`MultInt/ptorpols2`(eq,intnvars,degR,tdegR,strlst):

 if gu=0 then RETURN(NULL) fi:

fi:

ope:=gu[1]:

Rlst:=[]:

for i from 1 to nops(intnvars) do

 Rlst:=[op(Rlst), op(i,op(2,gu))/op(i,polylst)]:

od:

ope:=normal(ope):

mekh:=denom(ope):

ope:=numer(ope):

ope:=`MultInt/pashet2`(ope,mainvars,intnvars,fnam):

for i from 1 to nops(intnvars) do

 Rlst:=subsop(i=factor(Rlst[i]*mekh/p), Rlst):

od:

if nops(intnvars) >1 and  `MultInt/sym`(p*f,intnvars) then

 SR := `MultInt/symR`(Rlst,intnvars):
 
 printf(`%s%s%.2f%s\n`,` Cheers! for the success. `,  `  CPU Time : `,
            
             time() -st, ` seconds.`): print(``):
          
  
fvars:=op(intnvars):
mvars:=op(mainvars):
rhside:=0: 
for i from 1 to nops(intnvars) do 
deltintnvar:=op(subsop(i=NULL,intnvars)):
rhside:=rhside + Diff(_R(intnvars[i], deltintnvar)*fnam(mvars, fvars),intnvars[i]):
od:             
   
RETURN([ope=rhside, _R(fvars) = SR]):
   
#RETURN([ope=sum('Diff(_R(intnvars[i],op(subsop(i=NULL,intnvars)))*
 #            fnam(op(mainvars),op(intnvars)),intnvars[i])','i'=1..nops(intnvars)),
  #           _R(op(intnvars)) = SR]):

fi:
 
printf(`%s%s%.2f%s\n`,` Cheers! for the success. `,  ` CPU Time : `,
            
             time() - st, ` seconds.`): print(``):

fvars:=op(intnvars):
mvars:=op(mainvars):
rhside:=0: 
for i from 1 to nops(intnvars) do 
deltintnvar:=op(subsop(i=NULL,intnvars)):
rhside:=rhside + Diff(Rlst[i]*fnam(mvars, fvars),intnvars[i]):
od:             
   
RETURN(ope=rhside):

#RETURN(ope=sum('Diff(Rlst[i]*fnam(op(mainvars),op(intnvars)),intnvars[i])',
 #     'i'=1..nops(intnvars))):

end:

#---------------------------------------------------------------------------
`MultInt/findrec301`:=proc(p,f,mainvars,intnvars,axvars,strlst,fnam,st)

local i,j,den,den1,den11,den12,den2,factrs,factrs2,factr,indx,indxlst,result,
    
      polylst,maxd,maxord,dfbyf,denlst:

result:=NULL: den:=`MultInt/extrctdenom`(f,intnvars):

maxd:=0: factrs:=op(2,factors(den)):

for i from 1 to nops(factrs) do
 
 factr:=collect(op(1,op(i,factrs))^op(2,op(i,factrs)),intnvars,'distributed'):

 maxd:=max(maxd,degree(factr,{op(intnvars)})):
 
od:


if nops(strlst)<>0 then

 maxord:=max(op(map(proc(x) convert(x,`+`) end,strlst))):

else
 
 maxord:=0:
  
fi:

if maxord=0  then maxd:=0 fi:

den1:=select(type,den,{name ^ nonnegint, name ^ 1}):

if  result=NULL and den1<>NULL and has(den1,intnvars) then

 for i from 0 to maxord + maxd while result=NULL do

  den11:=factor(den1^i):
  
  polylst:=[seq(den11,j=1..nops(intnvars))]:

 
result:=`MultInt/findrec30`(p,f,mainvars,intnvars,axvars,strlst,polylst,fnam,st):

 od:

fi:

if result=NULL and den1<>NULL and has(den1,intnvars) then

 denlst:=[]:

 for i from 1 to nops(intnvars) do

  dfbyf:=diff(simplify(log(f)),intnvars[i]):  

  dfbyf:=normal(dfbyf):

  denlst:=[op(denlst),denom(dfbyf)]:

 od:

 polylst:=[seq(normal(op(i,denlst)/gcd(den1,op(i,denlst))),i=1..nops(denlst))]:
   
result:=`MultInt/findrec30`(p,f,mainvars,intnvars,axvars,strlst,polylst,fnam,st):

fi: 


if result=NULL and den1<>NULL and has(den1,intnvars) then

  den2:=factor(den/den1):

  factrs:=op(2,factors(den2)): factrs2:=[]:
 
  for i from 1 to nops(factrs) while result=NULL do

   factr:=op(i,factrs):
 
   factr:=op(1,factr)^(op(2,factr)): factr:=factor(factr):

   if has(factr,intnvars) then

    factrs2:=[op(factrs2),factr]:
 
    polylst:=[seq(factor(factr),j=1..nops(intnvars))]:
    
 result:=`MultInt/findrec30`(p,f,mainvars,intnvars,axvars,strlst,polylst,fnam,st):

     fi:

  od:

  if nops(factrs2)>1 and result=NULL then
  
    with(combinat):
   
   indxlst:=choose(nops(factrs2),2):
   
   for i from 1 to nops(indxlst) while result=NULL do
   
    indx:=op(i,indxlst):
    
    factr:= op(indx[1],factrs2)*op(indx[2],factrs2):
    
    factr:=factor(factr):
    
    polylst:=[seq(factr,j=1..nops(intnvars))]:

   result:=`MultInt/findrec30`(p,f,mainvars,intnvars,axvars,strlst,polylst,fnam,st):

   od:
   
  fi:
 
  for i from 0 to maxord + maxd while result=NULL do

    den12:=factor(den1^i*den2):
    
    polylst:=[seq(den12,j=1..nops(intnvars))]:

   result:=`MultInt/findrec30`(p,f,mainvars,intnvars,axvars,strlst,polylst,fnam,st):

   od:
 
fi:
 
if result=NULL and den1<>NULL and has(den1,intnvars) then

  for i from 2 to maxord + maxd while result=NULL do
  
   den12:=factor(den1*den2^i):
  
   polylst:=[seq(den12,j=1..nops(intnvars))]:

  result:=`MultInt/findrec30`(p,f,mainvars,intnvars,axvars,strlst,polylst,fnam,st):

  od: 

fi:
 
if result=NULL then
 
   for i from 0 to maxord + maxd while result=NULL do

    den1:=factor(den^i):
 
    polylst:=[seq(den1,j=1..nops(intnvars))]:

   result:=`MultInt/findrec30`(p,f,mainvars,intnvars,axvars,strlst,polylst,fnam,st):

   od:
   
fi:

if result=NULL then

 print(`  could not find a non-zero WZ eq with the given structure list`):            
fi:

result:
 
end:  

#--------------------------------------------------------------------------

`MultInt/findrec3`:=proc()

local ansatz,args1,nargs1,args2,i,p,f,strlst,polylst,t1,t2,shvars,set1,setdf,st:

if nargs < 6 then ERROR(` Wrong number of arguments`) fi:

#defining the correct args types

t1:=[algebraic,list(name),list(name),list(name),list(name),

    list({1,name,`*`,`^`}), list(polynom),name,numeric]:

t2:=[algebraic,list(name),list(name),list(name),list(name),
 
    list({1,name,`*`,`^`}),name,numeric]:

args1:=[args]: 


if nops(args1[2])<>nops(args1[3]) or 

   nops({op(args1[2])} intersect {op(args1[3])})<>0 then

  ERROR(cat(` the second and third arguments must have the same length\n`,

          ` and no common variable name. `)):
 
fi:

ansatz:=args1[6]: shvars:=args1[3]:

set1:={seq(type(ansatz[i],polynom(integer,shvars)),i=1..nops(ansatz))}:

if set1<>{true} then ERROR(cat(` invalid list `, convert(ansatz,string))) fi:

strlst:=`MultInt/structurelst`(ansatz,shvars):

nargs1:=nops(args1): 

args2:=subsop(1=NULL,nargs1-1=NULL,nargs1=NULL,args1):

if type(args1,t1) then

 polylst:=`MultInt/checkpoly`(args2[nops(args2)],nops(args1[4])):

 p:=`MultInt/extrctP`(args1[1],args1[2],args1[4],args1[5]):

 f:=p[2]: p:=p[1]: 

 args2:=subsop(2=NULL,nops(args2)=NULL,nops(args2)-1=NULL,args2):

 `MultInt/findrec30`(p,f,op(args2),strlst,polylst,args1[8],args1[9]):
 
elif type(args1,t2) then

 p:=`MultInt/extrctP`(args1[1],args1[2],args1[4],args1[5]):

 f:=p[2]: p:=p[1]:

 args2:=subsop(2=NULL,nops(args2)=NULL,args2):

 `MultInt/findrec301`(p,f,op(args2),strlst,args1[7],args1[8]):

fi:
 
end:


#-----------------------------------------------------------------------------
# `MultInt/checkrec`: takes a recurrence-differential (WZ-equation)
#                  returns True if f satisfies this WZ-equation.

`MultInt/checkrecdiff`:=proc(rec,f)

local t1,t2:

t1:=[`=`,algebraic]: t2:=[list(`=`),algebraic]:

if type([args],t1) then
 
  RETURN(`MultInt/checkrec1`(rec,f)):

fi:

if type([args],t2) then
 
 RETURN( `MultInt/checkrec2`(rec,f)):

fi:

ERROR(cat(convert([args],string),`\nmust be of type:\n`,

`[rec_diff equation, alg. expression]\n or `,

`[[rec_diff equation,_R[...]=rational expr.], alg. expression]`)):

end:

MultInt[checkrecdiff]:=proc()

`MultInt/checkrecdiff`(args):

end:

#------------------------------------------------------------------------------
`MultInt/checkrec1`:=proc(rec,f)

local expr,indets1,indets2,i,convars,ope,Difflst,Dpart,
      recprt,Dprt,s,fnam,f1,flst,flsti:

indets1:=indets(rec): convars:=[]:

indets1:=select(type, indets1, function):

if nops(indets1)=0 then 

 ERROR(` invalid rec_diff equation`)

fi:

Difflst:=convert(select(has, indets1, Diff),list):

if nops(Difflst)=0 then 

 ERROR(` invalid rec_diff equation`)

fi:

for i from 1 to nops(Difflst) do

convars:=[op(convars),op(2,Difflst[i])]:

od:

convars:=[op(convars)]:
 
indets2:=remove(has,indets1,Diff):

if nops(indets2)=0 then 

 ERROR(` invalid rec_diff equation`)

fi:

indets2:=convert(indets2,list):

flst:=[op(map(proc(x) op(0,x) end,indets2))]:

if nops(flst) = 1 and 

 {op(op(indets2))} intersect {op(convars)} = {op(convars)}  then

 fnam:=op(0,op(indets2)):

else

 fnam:=0:

 for i from 1 to nops(flst) while fnam=0 do

  flsti:= select(has,indets2,flst[i]):

  f1:=op(1,flsti):

  if {op(f1)} intersect {op(convars)} = {op(convars)} then
  
   fnam:=op(0,f1):

  fi:
 
 od:

 if fnam=0 then

  ERROR(` invalid rec_diff equation`)

 fi:

fi:

ope:=lhs(rec): Dpart:=rhs(rec):  

if has(ope,Diff) then
  
 expr:=select(has,ope,Diff) :
   
 if has(expr,Diff) then 
   
  ope:=ope - expr: 
    
  Dpart:=Dpart+expr:
    
 fi:
   
fi:

if has(Dpart,fnam) then
  
  expr:=select(has,Dpart,fnam) :
   
  if has(expr,fnam) and not has(expr,Diff) then 
   
   ope:=ope + expr: 
    
   Dpart:=Dpart - expr:
    
  fi:
   
fi:  

recprt:=`MultInt/evalrecprt`(ope,f,fnam):

Dprt:=`MultInt/evalDprt1`(Dpart,f,fnam):

s:=simplify(recprt-Dprt):

if s<>0 then
 
 false
 
else

 true

fi:

end:

#------------------------------------------------------------------------------
`MultInt/evalrecprt`:=proc(ope,f,fnam)

local summands, summand,i,j,n,v1,var, v2, sh, f1, indx,ope1:

ope1:=0:

if type(ope,`+`) then

  summands:=[op(ope)]:
  
  for i from 1 to nops(summands) do
  
   summand:=summands[i]:
   
   sh:=convert(select(has,indets(summand),fnam),list):

   if not (nops(sh)=1 and op(0,op(sh))=fnam) then 
   
     ERROR(` invalid rec_diff equation. `)

   fi:

   sh:=op(sh):

   summand:=summand/sh:
   
   indx:=[op(sh)]:
   
   f1:=f:
   
   for j from 1 to nops(indx) do
       
    v1:=indx[j]:
   
    if type(v1,`+`) then
   
     v2:={op(v1)}:
    
     n:=select(type,v2,numeric): var:=select(type,v2,name):
    
     if nops(n)=1 and nops(var)=1  then
     
       n:=op(n): var:=op(var):
       
       f1:=subs(var=var +n,f1)      
       
     elif nops(n)>1 or nops(var)>1 then
    
       ERROR(` Invalid rec_diff equation `):
     
     fi:
     
    fi:
  
   od:
  
   f1:=simplify(summand*f1/f):
  
   ope1:=ope1 + f1:
 
  od:
 
 else
     
   sh:=convert(select(has,indets(ope),fnam),list):

   if not (nops(sh)=1 and op(0,op(sh))=fnam) then 
   
     ERROR(` invalid rec_diff equation. `)

   fi:

   sh:=op(sh):

   summand:=ope/sh:
   
   indx:=[op(sh)]:
   
   f1:=f:
   
   for j from 1 to nops(indx) do
   
    v1:=indx[j]:
   
    if type(v1,`+`) then
   
     v2:={op(v1)}:
    
     n:=select(type,v2,numeric): var:=select(type,v2,name):
    
     if nops(n)=1 and nops(var)=1  then
     
       n:=op(n): var:=op(var):
       
       f1:=subs(var=var +n,f1)      
       
     elif nops(n)>1 or nops(var)>1 then
    
       ERROR(` Invalid rec_diff equation `):
     
     fi:
     
    fi:
  
   od:
  
   f1:=simplify(summand*f1/f):
  
   ope1:=ope1 + f1:
 
 fi:

RETURN(ope1):
  
end:



#-----------------------------------------------------------------------------

`MultInt/evalDprt1`:=proc(Dpart,f,fnam)

local summands, summand,i,var, R, RF, Fprt, s:

s:=0:

if type(Dpart,`+`) then

  summands:=[op(Dpart)]:
  
  for i from 1 to nops(summands) do
  
   summand:=summands[i]:
   
   if not (type(summand,function) and op(0,summand)=Diff) then
   
     ERROR(` invalid rec_diff equation`):
     
   fi:
   
   var:=op(2,summand): RF:=op(1,summand):
   
   if RF<>0 then

    Fprt:=select(has,RF,fnam):
   
    R:=RF/Fprt:
   
    s:=s + factor(diff(simplify(log(R*f)),var)*R):
    
   fi:
    
  od:     
  
elif type(Dpart,function) and op(0,Dpart)=Diff then
    
   var:=op(2,Dpart): RF:=op(1,Dpart):
   
   Fprt:=select(has,RF,fnam):
   
   R:=RF/Fprt:
   
   s:=s + factor(diff(simplify(log(R*f)),var)*R):
   
else

  ERROR(` Invalid rec_diff equation`)
  
fi:

RETURN(s):
  
end:

#------------------------------------------------------------------------------

`MultInt/checkrec2`:=proc(rec,f)

local expr,indets1,i,indets2,Difflst,fnam,f1,flst,flsti,convars,
      
      intnvars,receq,ope,Dpart,recprt,Dprt,s,RC:

receq:=rec[1]: RC:=rhs(rec[2]): 

if op(0,lhs(rec[2]))=_R and 

  map(proc(x) type(x,name) end,{op(lhs(rec[2]))})={true} then

 intnvars:=[op(lhs(rec[2]))]:

else

 ERROR(` invalid recurrence-diff. equation `,rec):
 
fi:


indets1:=indets(receq): convars:=[]:

indets1:=select(type, indets1, function):

if nops(indets1)=0 then 

 ERROR(` invalid rec_diff equation`)

fi:

Difflst:=convert(select(has, indets1, Diff),list):

if nops(Difflst)=0 then 

 ERROR(` invalid rec_diff equation `,rec)

fi:

for i from 1 to nops(Difflst) do

convars:=[op(convars),op(2,Difflst[i])]:

od:

convars:=[op(convars)]:
 
indets2:=remove(has,indets1,{Diff,_R}):

if nops(indets2)=0 then 

 ERROR(` invalid recurrence-diffferential equation`)

fi:

indets2:=convert(indets2,list):

flst:=[op(map(proc(x) op(0,x) end,indets2))]:

if nops(flst) = 1 and 

 {op(op(indets2))} intersect {op(convars)} = {op(convars)}  then

 fnam:=op(0,op(indets2)):

else

 fnam:=0:

 for i from 1 to nops(flst) while fnam=0 do

  flsti:= select(has,indets2,flst[i]):

  f1:=op(1,flsti):

  if {op(f1)} intersect {op(convars)} = {op(convars)} then
  
   fnam:=op(0,f1):

  fi:
 
 od:

 if fnam=0 then

  ERROR(` invalid recurrence-diff equation`)

 fi:

fi:

ope:=lhs(receq): Dpart:=rhs(receq):

if has(ope,Diff) then
  
 expr:=select(has,ope,Diff) :
   
 if has(expr,Diff) then 
   
  ope:=ope - expr: 
    
  Dpart:=Dpart+expr:
    
 fi:
   
fi:

if has(Dpart,fnam) then
  
  expr:=select(has,Dpart,fnam) :
   
  if has(expr,fnam) and not has(expr,Diff) then 
   
   ope:=ope + expr: 
    
   Dpart:=Dpart - expr:
    
  fi:
   
fi:  


recprt:=`MultInt/evalrecprt`(ope,f,fnam):

Dprt:=`MultInt/evalDprt2`(Dpart,RC,intnvars,f,fnam):

s:=simplify(recprt-Dprt):

if s<>0 then
 
 false
 
else

 true

fi:

end:

#-----------------------------------------------------------------------------
`MultInt/evalDprt2`:=proc(Dpart,RC,intnvars,f,fnam)

local summands, summand,i,j,var, R, Rvars, RF, Fprt, s:

s:=0:

if type(Dpart,`+`) then

  summands:=[op(Dpart)]:
  
  for i from 1 to nops(summands) do
  
   summand:=summands[i]:
   
   if not (type(summand,function) and op(0,summand)=Diff) then
   
     ERROR(` invalid rec_diff equation`):
     
   fi:
   
   var:=op(2,summand): RF:=op(1,summand):
   
   if RF<>0 then

     Fprt:=select(has,RF,fnam):
   
     R:=RF/Fprt:
   
    if op(0,R)=_R then
      
     Rvars:=[op(R)]:
    
     R:=subs({seq(intnvars[j]=Rvars[j],j=1..nops(Rvars))},RC):   
   
     s:=s + factor(diff(simplify(log(R*f)),var)*R):
    
    else
    
     ERROR(` Invalid recurrence-differential equation`):
    
    fi:
    
   fi:
   
  od:     
  
elif type(Dpart,function) and op(0,Dpart)=Diff then
    
   var:=op(2,Dpart): RF:=op(1,Dpart):
   
   Fprt:=select(has,RF,fnam):
   
   R:=RF/Fprt:
   
   s:=s + factor(diff(simplify(log(R*f)),var)*R):
   
else

  ERROR(` Invalid rec_diff equation`)
  
fi:

RETURN(s):
  
end:

`MultInt/rf`:=proc(x,m)
 
GAMMA(x+m)/GAMMA(x):
 
end:

#-------------------------------------------------------------------------------
`MultInt/extrctnCk`:=proc(lst,k)

local aCb, nCk, a, b, i, lst1,pos:

aCb:=op(1,lst): a:=op(1,aCb): b:=op(2,aCb):

if degree(b,k)=1 then pos:=1 else pos:=0 fi:

for i from 1 to nops(lst) do

 aCb:=op(i,lst):

 a:=op(1,aCb): b:=op(2,aCb):
 
 if not has(a,k) and degree(b,k)=1 then
  
   if pos=0 then pos:=i: fi:
   
   if nops([coeffs(b,k)])<nops([coeffs(op(2,lst[pos]),k)])  then
    
       pos:=i:

   fi:
   
 fi:

od:

if pos=0 then RETURN([],lst) fi:

nCk:=op(pos,lst):

lst1:= subsop(pos=NULL,lst):

nCk,lst1:

end:

#------------------------------------------------------------------------------
`MultInt/powerofBC`:=proc(BC,var)

local b,expon,result,i:

if type(BC,`^`) then

 b:=op(1,BC): expon:=op(2,BC):   if not (type(b,function) and op(0,b)=binomial) then    ERROR(` Invalid expression: `, b):   fi:     b:=[op(1,b),op(2,b)]:   if type(b,[polynom(anything,var),polynom(anything,var)]) then

   if type(expon,posint) then

     result:=[seq(b,i=1..expon)]:

   else

    ERROR(` Binomial coeff. with negative power.`):

   fi:
   
 else
 
   ERROR(` Invalid Binomial coeff: `,BC):
 
 fi:
 
else
 
 ERROR(` Invalid Binomial coeff: `,BC):

fi:

result:

end:

#-----------------------------------------------------------------------------
`MultInt/makbclst`:=proc(BC,var)

local i,t,lst,lstBC,b:


t:={`^`,function,`*`}:

if not type(BC,t) then

 ERROR(` Invalid expression `,BC):

fi:

lstBC:=[]:

if type(BC,`^`) then

  lstBC:=[op(lstBC), op(`MultInt/powerofBC`(BC,var))]:

elif type(BC,`*`) then

 lst:=[op(BC)]:
 
 for i from 1 to nops(lst) do
 
  if type(op(i,lst),`^`) then
   
   lstBC:=[op(lstBC), op(`MultInt/powerofBC`(op(i,lst),var))]:
   
  elif type(op(i,lst),function) and op(0, op(i,lst))=binomial then
       b:=[op(1,op(i,lst)),op(2,op(i,lst))]:       if type(b,[polynom(anything,var),polynom(anything,var)]) then
     lstBC:=[op(lstBC),b]:     else   
    ERROR(` Invalid expression in `, binomial(op(1,b),op(2,b))):     fi:    
  else
  
   ERROR(` Invalid binomial coeff expression.`):
   
  fi:
   
 od:
 
elif type(BC,function) and op(0,BC)=binomial  then   b:=[op(1,BC),op(2,BC)]:   if type(b,[polynom(anything,var),polynom(anything,var)]) then     lstBC:=[op(lstBC),b]:   else  
   ERROR(` Invalid expression `,binomial(op(1,b),op(2,b))):   fi: 
fi:

lstBC:

end:
 
#---------------------------------------------------------------------------

`MultInt/splitlst`:=proc(lst,k)

local i,j,lst1, lsti, nCk, indx, klst, klst2, klst3, k1,x,pos:

lst1:=lst: klst:=[]:

for i from 1 to nops(k) do
  
  lsti:=[]:
  
  indx:=op(i,k): 
  
  for j from 1 to nops(lst1) do
  
    nCk:=op(j,lst1):
    
    if not has(op(1,nCk),{indx}) and degree(op(2,nCk),indx)=1 then
    
     lsti:=[op(lsti),nCk]:
     
    fi:
    
  od:   
  

  lst1:={op(lst1)} minus {op(lsti)}:
     
  if nops(lsti)<>0 then
  
    klst:=[op(klst),[indx,lsti]]:
 
  fi:
   
od:

klst2:=[]: klst3:=[]:
  
if nops(klst) <> 0 then

 for i from 1 to nops(klst) do
 
  lsti:=op(i,klst):
 
  k1:=op(1,lsti):  lsti:=op(2,lsti):
  
  nCk:=`MultInt/extrctnCk`(lsti,k1):

  if nops(nCk[1])<>0 then
  
   klst2:=[op(klst2),nCk[1]]:
   
   klst3:=[op(klst3),[k1,nCk[1]]]:
      
  fi:
  
 od:
 
fi:

lst1:=lst:

for i from 1 to nops(klst2) do

 x:=op(i,klst2): pos:='pos':
 
 if member(x,lst1,`pos`) then
  
  lst1:=subsop(pos=NULL,lst1):
 
 fi:

od:

klst3,lst1:
  
end:

#----------------------------------------------------------
`MultInt/sumtointn`:=proc()

local t1,t2:

if nargs<4 then

 ERROR(` wrong number of arguments `):
 
fi:


t1:=[{`*`,`^`,function},algebraic,name, [algebraic,algebraic]]:

t2:=[{`*`,`^`,function},algebraic, list(name), list([algebraic,algebraic])]:

if type([args], t1) then


 if `MultInt/notdistvars`(op(3,[args])) then

   ERROR(cat(` not distinct variables: `, convert(op(3,[args]),string))):
  
 fi:

 RETURN(`MultInt/ssumtointn`(args)):

elif type([args], t2) then

 if `MultInt/notdistvars`(op(3,[args])) then

  ERROR(cat(` not distinct variables: `, convert(op(3,[args]),string))):
   
 fi:

 if nops(op(3,[args]))<>nops(op(4,[args])) then

  ERROR(cat(` there should be 1-1 correspondence between `,
   
   convert(op(3,[args]),string), ` and `,

   convert(op(4,[args]),string))):

 fi:    

 
 RETURN(`MultInt/msumtointn`(args)):       

else
 

 ERROR(cat(` expecting argument of type:\n`,
 
  `(finite product of) binomial coeff(s), x, summation var, `, 
  
   ` summation bounds `, `\n or \n`, 

   ` product of binomial coeffs., x, list of summation vars, `,
  
  ` list of intervals of sum bounds`)):
  
fi:

end:

MultInt[sumtointn]:=proc()

`MultInt/sumtointn`(args):

end:


#-----------------------------------------------------------------------------    
  
`MultInt/ssumtointn`:=proc(bc,x,k,interval)

local t,lc,ct,nCk,expon,lst,aCb,a,b,prd,i,k1,n1,result,v:

lst:=`MultInt/makbclst`(bc,k):

nCk:=`MultInt/extrctnCk`(lst,k): 

lst:=nCk[2]: nCk:=nCk[1]: 

prd:=1: 

for i from 1 to nops(lst) do  
  
  aCb:=op(i,lst):   
  
  a:=op(1,aCb): b:=op(2,aCb):
   
  prd:=prd*(1+1/(cat(z_,i)))^a*(cat(z_,i))^b:
  
od:

prd:=prd*x:

a:=op(1,interval):

b:=op(2,interval):

if nops(nCk)<>0 then

 n1:=op(1,nCk): k1:=op(2,nCk):
 
 lc:=lcoeff(k1,k): ct:=coeff(k1,k,0):
 
 prd:=eval(subs(k=(k-ct)/lc,prd)):
  
 n1:=eval(subs(k=(k-ct)/lc,n1)):
 
 k1:=eval(subs(k=(k-ct)/lc,k1)):

 result:= sum(binomial(n1,k1)*prd,k=0..infinity):

else

 result:=sum(prd,k=a..b):   
 
fi:

v:=[seq(cat(z_,i),i=1..nops(lst))]:

if has(result,binomial) then

 result:=convert(result,factorial):

fi:

result:=normal(result):

if nops(v)<>0 then
 
 RETURN(result, v):
 
else

 RETURN(result,[cat(z_,1)]):

fi:
 
end:


#-----------------------------------------------------------------------------

`MultInt/newsum`:=proc(prd,lst,k,interval,v)

local i,nCk,pos,lst1,nv, prd1,lst2,lc,ct,k1,k2,newk,poslst,poslst2,interval1,v1:

prd1:=prd: lst1:=lst:  lst2:=[]: newk:=[]: poslst:=[]: v1:=v:

nv:=nops(v):
  
for i from 1 to nops(lst) do
  
 nCk:=op(i,lst1):
 
 k1:=op(1,nCk): nCk:= op(2,nCk):
 
 pos:='pos': newk:=[op(newk),k1]:
 
 if member(k1,k,pos) then
 
  poslst:=[op(poslst),pos]:
  
 fi:
 
 lst2:=[op(lst2),nCk]:
 
 k2:=op(2,nCk):
  
 lc:=lcoeff(k2,k1): ct:=coeff(k2,k1,0):
 
 if lc<>1 or ct<>0 then
  
  prd1:=eval(subs(k1=(k1-ct)/lc,prd1)):
   
  lst1:=eval(subs(k1=(k1-ct)/lc,lst1)):
  
  lst2:=eval(subs(k1=(k1-ct)/lc,lst2)):
  
 fi:  
 
od:

poslst2:=[]:

for i from 1 to nops(lst2) do

 nCk:=op(i,lst2): 
  
 if has(op(1,nCk),op(2,nCk)) then
 
  nv:=nv+1:
 
  prd1:=prd1*(1+1/(cat(z_,nv)))^op(1,nCk)*(cat(z_,nv))^op(2,nCk):
  
  v1:=[op(v1),cat(z_,nv)]:
  
  poslst2:=[op(poslst2),i]:
  
 fi:  
    
od: 

k1:=subsop(seq(op(i,poslst)=NULL,i=1..nops(poslst)),k):

interval1:=subsop(seq(op(i,poslst)=NULL,i=1..nops(poslst)),interval):
   
if nops(poslst2)<>0 then

 lst2:= subsop(seq(op(i,poslst2)=NULL,i=1..nops(poslst2)),lst2):
 
 k1:=[op(k1),seq(op(op(i,poslst2),newk),i=1..nops(poslst2))]:
 
 interval1:=[op(interval1),seq([0,infinity],i=1..nops(poslst2))]:
 
 newk:= subsop(seq(op(i,poslst2)=NULL,i=1..nops(poslst2)),newk):
 
fi: 
   

prd1,lst2,newk, k1, interval1,v1:

end:

#-----------------------------------------------------------------------------

`MultInt/orderindx`:=proc(blst,vlst)
 
local i, vlst1,blst1,v,found:

vlst1:=vlst: blst1:=blst: v:=[]: found:=false:
 
for i from 1 to nops(blst1) do
 
 v:=[op(v),op(1,op(i,blst1))]:

od:
  
for i from 1 to nops(vlst) while found=false do
 
 if not has(v, vlst[i]) then
  
  found:=true:
    
  blst1:=subsop(1=op(i,blst1),i=op(1,blst1),blst1):
   
  vlst1:=subsop(1=op(i,vlst1),i=op(1,vlst1),vlst1):
   
 fi:
   
od:
 
blst1, vlst1:

end:
  
#-------------------------------------------------------------------------------
`MultInt/msumtointn`:=proc(bc,x,k,interval)

local nCk, lst,lst1, aCb, a, b,t, prd,lstk,i,k1,n1,s,indx,newk1,newk2,

     interval1,expr, v:

lst:=`MultInt/makbclst`(bc,k):

lst1:=`MultInt/splitlst`(lst,k):

lstk:=lst1[1]: lst1:=lst1[2]: 

prd:=1:
 
for i from 1 to nops(lst1) do

 aCb:=op(i,lst1): a:=op(1,aCb): b:=op(2,aCb):
 
 prd:=prd*(1+1/(cat(z_,i)))^a*(cat(z_,i))^b:
   
od:

v:=[seq(cat(z_,i),i=1..nops(lst1))]:

prd:=prd*x: 

if nops(lstk) <> 0 then

 expr:=`MultInt/newsum`(prd,lstk,k,interval,v):
 
 prd:=expr[1]: lstk:=expr[2]: newk1:=expr[3]: 
 
 newk2:=expr[4]: interval1:=expr[5]: v:=expr[6]:
 
 lstk:= `MultInt/orderindx`(lstk,newk1):
 
 newk1:=lstk[2]: lstk:=lstk[1]: 
 
 nCk:=op(1,lstk):

 n1:=op(1,nCk): k1:=op(2,nCk): indx:=op(1,newk1): 
 
 s:=sum(binomial(n1,k1)*prd,indx=0..infinity):
 
 newk1:=subsop(1=NULL,newk1): lstk:=subsop(1=NULL,lstk):
  
 while nops(lstk)<>0 do
 
  lstk:= `MultInt/orderindx`(lstk,newk1):
 
  newk1:=lstk[2]: lstk:=lstk[1]: 
   
  nCk:=op(1,lstk): 

  n1:=op(1,nCk): k1:=op(2,nCk): indx:=op(1,newk1):
  
  s:=sum(binomial(n1,k1)*s,indx=0..infinity):
  
  newk1:=subsop(1=NULL,newk1): lstk:=subsop(1=NULL,lstk):
   
 od:

 for i from 1 to nops(newk2) do
   
   a:=op(1,interval1[i]):

   b:=op(2,interval[i]):
  
   k1:=op(i,newk2):

   s:=sum(s,k1=a..b):   
   
 od:
  
else

 interval1:=op(1,interval):
  
 a:=op(1,interval1):

 b:=op(2,interval1):

 k1:=op(1,k):
 
 s:=sum(prd,k1=a..b):   
 
 for i from 2 to nops(interval) do

  interval1:=op(i,interval):
  
  a:=op(1,interval1):

  b:=op(2,interval1):
  
  k1:=op(i,k):

  s:=sum(s,k1=a..b):   
  
 od:
 
fi:

if has(s,binomial) then

  s:=convert(s,factorial):
  
fi:

s:=normal(s):
 
if nops(v)<>0 then
 
 RETURN(s, v):
 
else

 RETURN(s,[cat(z_,1)]):

fi:
 
end:

#-------------------------------------------------
`MultInt/sumtorecdiff`:=proc()

local args1,t1,t2,t3,t4:

if nargs < 5 then ERROR(` Wrong number of arguments. `) fi:

#defining the correct args types

args1:=[args]:

t1:=[{`*`,`^`,function}, algebraic, name,
     
     name, [algebraic,algebraic], nonnegint]:

t2:=[{`*`,`^`,function}, algebraic, name, 
     
     name, [algebraic,algebraic]]:


 
t3:=[{`*`,`^`,function}, algebraic, list(name),
     
     name, list([algebraic,algebraic]), nonnegint]:

t4:=[{`*`,`^`,function}, algebraic, list(name), 
     
     name, list([algebraic,algebraic])]:


if type(args1,{t1,t2}) then

 if `MultInt/notdistvars`([args1[3],args1[4]]) then

  ERROR(cat(` not distinct variables: `,
  
   convert([args1[3],args1[4]],string))):
  
 fi:

 if nops([op(3,args1)])<>nops([op(5,args1)]) then

  ERROR(cat(` there should be 1-1 correspondence between `,
   
    op(3,args1), ` and `, convert(op(5,args1),string))):

 fi:    
       
 if not type([op(3,args1),op(4,args1)],[name,name]) then

  ERROR(cat(` invalid argument:\n`, 
 
     convert([op(3,args1),op(4,args1)],string))) :

 fi:

 RETURN(`MultInt/ssum`(args)):

elif type(args1,{t3,t4}) then

 if `MultInt/notdistvars`([op(args1[3]),args1[4]]) then

  ERROR(cat(` not distinct variables: `, convert([op(args1[3]),args1[4]],string))):
  
 fi:

 if nops(op(3,args1))<>nops(op(5,args1)) then

 ERROR(cat(` there should be 1-1 correspondence between `,
   
    convert(op(3,args1),string), ` and `, convert(op(5,args1),string))):

 fi:    
       
 if not type([op(3,args1),op(4,args1)],[list(name),name]) then

  ERROR(cat(` invalid argument:\n`,

  convert([op(3,args1),op(4,args1)],string))) :

 fi:

 RETURN(`MultInt/msum`(args)):

else

 ERROR(cat(` Invalid input, expecting an arg type:\n `, 
 
 `(finite product of) binomial coeffs, x , sum var, recurrence var`,
 
 ` sum bound `,
 
  `\nor\n `, `(finite product of) binomial coeffs, x , sum var, recurrence var`,

        ` sum  bound`, ` recurrence order`, `\nor\n `,

   `(finite product of) binomial coeffs, x, list of sum vars, recurrence var`,  

  ` list of sum bounds `,
 
  `\nor\n `, ` (finite product of) binomial coeffs, x , list of sum vars`,
  
  ` recurrence var, list of sum bound`, ` recurrence order`)):

fi:

end:

MultInt[sumtorecdiff] := proc()

`MultInt/sumtorecdiff`(args):

end:


#--------------------------------------------------------------------------

`MultInt/ssum`:=proc()

local args1,args2,i,t1,t2,axvars,Fz,vars, F, result:

args1:=[args]:

t1:=[{`*`,`^`,function}, algebraic, name,
     
     name, [algebraic,algebraic], nonnegint]:

t2:=[{`*`,`^`,function}, algebraic, name, 
     
     name, [algebraic,algebraic]]:


if type(args1,t1) then   args2:=seq(op(i,args1),i=1..3):

 Fz:=`MultInt/ssumtointn`(args2,args1[5]):
 
 F:=Fz[1]: vars:=Fz[2]: 
 
 if has(F,hypergeom) then
  
  ERROR(cat(`MultInt, had difficulty in converting sum to integral`,
   
            `check your summation index bound`)):
 fi:
  
F:=normal(F/convert(vars,`*`)/(2*Pi*I)^nops(vars)):
  
result:=`MultInt/hypertorecdiff`(F,op(4,args1),_F(op(vars)),op(6,args1)):
 
if result<>0 then 
 
  if type(result,list) then
   
   result:=[op(result), _F(args1[4], op(vars)) = F ]:
   
  else
   
   result:=[result,_F(args1[4], op(vars)) = F]:
   
  fi:

fi:
 
 result:
  
elif type(args1,t2) then

 args2:=seq(op(i,args1),i=1..3):

 Fz:=`MultInt/ssumtointn`(args2,args1[5]):
 
 F:=Fz[1]: vars:=Fz[2]: 
 
 if has(F,hypergeom) then
  
   ERROR(cat(` MultInt, had difficulty in converting sum to integral`,
   
            `check your summation index bound.`)):
 fi:
 
 F:=normal(F/convert(vars,`*`)/(2*Pi*I)^nops(vars)):
 
 result:=`MultInt/hypertorecdiff`(F,op(4,args1),_F(op(vars))):
 
 if result<>0 then 
  
  if type(result,list) then
   
   result:=[op(result), _F(args1[4], op(vars)) = F ]:
   
  else
   
   result:=[result,_F(args1[4], op(vars)) = F]:
   
  fi:
  
 fi:
 
 result:

fi:

end:


#-------------------------------------------------------------------------------
`MultInt/msum`:=proc()

local args1,args2,i,t1,t2,axvars,Fz,vars, F,result:

args1:=[args]:
 
t1:=[{`*`,`^`,function}, algebraic, list(name),
     
     name, list([algebraic,algebraic]), nonnegint]:

t2:=[{`*`,`^`,function}, algebraic, list(name), 
     
     name, list([algebraic,algebraic])]:

if type(args1,t1) then

 args2:=seq(op(i,args1),i=1..3):

 Fz:=`MultInt/msumtointn`(args2,args1[5]):
  
 F:=Fz[1]: vars:=Fz[2]: 
 
 F:=normal(F/convert(vars,`*`)/(2*Pi*I)^nops(vars)):
 
 if has(F,hypergeom) then
  
  ERROR(cat(` MultInt, had difficulty in converting sum to integral`,
   
            ` check your summation index bound`)):
 fi:
 
 result:=`MultInt/hypertorecdiff`(F,op(4,args1),_F(op(vars)),op(6,args1)):
 
 if result<>0 then 
 
  if type(result,list) then
   
   result:=[op(result), _F(args1[4], op(vars)) = F ]:
   
  else
   
   result:=[result,_F(args1[4], op(vars)) = F]:
   
  fi:
  
 fi:
 
 result:
  
elif type(args1,t2) then

 args2:=seq(op(i,args1),i=1..3):

 Fz:=`MultInt/msumtointn`(args2,args1[5]):
 
 F:=Fz[1]: vars:=Fz[2]: 
 
 F:=normal(F/convert(vars,`*`)/(2*Pi*I)^nops(vars)):
 
 if has(F,hypergeom) then
  
  ERROR(cat(` MultInt, had difficulty in converting sum to integral`,
   
            ` check your summation index bound`)):
 fi:
 
 result:=`MultInt/hypertorecdiff`(F,op(4,args1),_F(op(vars))):
 
 if result<>0 then 
 
  if type(result,list) then
   
   result:=[op(result), _F(args1[4], op(vars)) = F ]:
   
  else
   
   result:=[result,_F(args1[4], op(vars)) = F]:
   
  fi:
  
 fi:
 
 result:

fi:

end:

mylatex:=proc(x,filename)
local `latex/special_names`:
readlib(latex):
`latex/special_names`['beginmath']:=
`\\begin{displaymath}`:
`latex/special_names`['endmath']:=
`\\end{displaymath}`:
latex(beginmath); latex(x,filename); latex(endmath);
end:


#_____________________________ HELP INFORMATION  ______________________________

print(` MultInt: a Maple Package for Multiple Integration `):
print(` of proper-hyperexponential functions by `):
print(` the continuous version of the multi-WZ Method.`):
print(``):
print(` Akalu Tefera, Grand Valley State University, Department of Mathematics.`):
print(``):
print(` Please report all bugs and comments to: teferaa@gvsu.edu`):
print(``):
print(` For a list of procedures, type: Help(); `): 
print(` `):
print(` For help with a specific procedure, type: Help(procedure name); `):
print(``):
print(` CAUTION: this version of MultInt is for Maple 6 and above `):
print(``):

if not member(`.`,[libname]) then

  libname:=libname,`.`:
  
fi:
  
#______________________________________________________________________________________

Help:=proc()
local rec:

if args=NULL then
print(`MultInt: a package for Mutliple Integration of proper-hyperexponential functions by the WZ method.`):
print(` `):
printf( "%s\n",` Calling  Sequences:    `): 
printf( "%s\n",`    MultInt[<function>](arguments) `):
printf( "%s\n", `   <function>(arguments)   `):

print(` `):

printf("%s\n",`Description: `):
printf("%s\n",`The functions available are:`):
print( `hypertorecdiff, checkrecdiff,esp, sumtointn,  sumtorecdiff , sym `):                 
printf("%s\n",`To get help for a specific function type Help(<function>) where <function> is one from the above list.`):   
printf("%s\n",`These functions are part of the MultInt package and so can be used in the form <function>(arguments) only after performing the command with(MultInt) or with(MultInt, <function>). The functions can also be accessed in the long form MultInt[<function>](arguments).`):
printf("%s\n",`Whenever there is a conflict between a function name in MultInt and another name used in the same session, use the form MultInt[<function>].`):
printf("%s\n",`This version of  MultInt is for Maple 6 and above.`):
fi:
    
if nops([args])=1 and op(1,[args])=`hypertorecdiff` then

print(`hypertorecdiff: finds a non-zero recurrence-differential(WZ) equation staisfied by a given proper-hyperexponential function.`):
print(` `):
printf("%s\n",`Calling Sequences:`):
printf("%s\n",`  hypertorecdiff(f, n, fnam(x1, x2, ..., xr), opt1, opt2,...)`):
printf("%s\n",`  hypertorecdiff(f, [n1, ..., nm],  fnam(x1, x2, ..., xr), opt1, opt2,...)`):
print(``):
printf("%s\n",`Paremeters: `):
printf("%s\n",`     f                - a proper-hyperexponential expressions`):
printf("%s\n",`     n, n1, ...,nm    - names, the recurrence variables`):
printf("%s\n",`     x1, x2,...,xr    - names, the integration variables`):
printf("%s\n",`     fnam             - a name, the function name`):
printf("%s\n",`     opt1, opt2, ...  - various desired options`):
print(` `):
printf("%s\n",`Description`):
printf("%s\n",`-hypertorecdiff(f, n, fnam(x1,x2,...,xr)) finds a non-zero recurrence-differential(WZ) equation of the form `):
printf("%s\n",`    p1*fnam(n,x1,..xr) + p2*fnam(n+1,x1,...,xr)+... = Diff(R1*fnam(n,x1,...,xr))+..., `):
printf("%s\n",`   where  p1p2,... are polynomials  independent (free) of  x1, ..., xr.`):
printf("%s\n",` In this case hypertorecdiff searches for a WZ-equation with recurrence part of order at most 6(the default value).`):
printf("%s\n",` - hypertorecdiff(f, [n1,...,nm], fnam(x1,x2,...,xr)) finds a non-zero multiple recurrence WZ-equation. `):
printf("%s\n",` The following options are supported:  ` ):
printf("%s\n",`- hypertorecdiff(f, n, fnam(x1,x2,...,xr), recurrence_order)`):
printf("%s\n",`    recurrence_order  -  non-negative integer`):
printf("%s\n",`    This option allows one  to input the required order of the recurrence part of the WZ-equation.`):
printf("%s\n",`    If n is a list of variables , then  recurrence_order is a list of non-negative integers.`):
printf("%s\n",` - hypertorecdiff(f, n, fnam(x1, x2, ..., xr), recurrence_order, denom_poly)`):
printf("%s\n",`    denom_poly -  list of  polynomials.`):
printf("%s\n",`    This option allows one to  guess and input the  denominators denom_poly  of the  r-tuple of rational functions(the certificates) of  the  required WZ-equation.`):
printf("%s\n",`    If n is a list of variables , then  recurrence_order is a list of non-negative integers.`):
printf("%s\n",` - hypertorecdiff(f, n, fnam(x1, ..., xr), denom_poly) `):
printf("%s\n",`    n can also be a list of recurrence variables. `):
printf("%s\n",` - hypertorecdiff(f, [n1,n2,...,nk], fnam(x1, ..., xr), [N1,...,Nk], ansatz)  `):
printf("%s\n",`    n1,....,nk - recurrence variables`):
printf("%s\n",`    N1,...Nk - shift variables`):
printf("%s\n",`    ansatz  - list of monomials in N1, ...,Nk`):
printf("%s\n",`    In this case hypertorecdiff searches for  a WZ-equation by using the given ansatz.`):
printf("%s\n",`- hypertorecdiff(f, [n1,n2,...,nk], fnam(x1, ..., xr), [N1,...,Nk], ansatz, denom_poly ) `):
printf("%s\n",`    This option allows the user to guess and input the denominators of the  r-tuple of rational functions(the certificates) of  the  required WZ-equation.`):
printf("%s\n",`This function is a part of the MultInt package, and so can be used in the form hypertorecdiff(args) only after performing the command with(MultInt) or  With(MultInt, hypertorecdiff). The function can also be accessed in the long form MultInt[hypertorecdiff](args).`):
printf("%s\n",`Whenever there is a conflict between  the function name hypertorecdiff and another name used in the same session, use the form MultInt['hypertprecdiff']`):

printf("%s\n",`Examples: `):
printf("%s\n",`> with(MultInt):`):
printf("%s\n",`> hypertorecdiff((1+x)^n*(1+y)^n/(1-x*y),n,f(x,y)); `):
 
print([(4*n+2)*f(n,x,y)+(-n-1)*f(n+1,x,y) = Diff(_R(x,y)*f(n,x,y),x)+Diff(_R(y,x)*f(n,x,y),y), _R(x,y) = -1/2*(1+x)*(x*y+2*x-3)]):

print(` `):
print(`The  following call  asks MultInt to find a non-zero WZ-equation of order 1.`):
printf("%s\n",`> hypertorecdiff((1+x)^n*(1+y)^m/(1-x*y)^k,m,f(x,y),1);`):
 print((n+m+2-k)*f(m,x,y)+(k-m-2)*f(m+1,x,y) = Diff((1+x)*f(m,x,y),x)+Diff(-y*(1+y)*f(m,x,y),y)):

print(` `):

print(` The  following calls  ask  MultInt to find a non-zero double recurrence  WZ-equation.`):

printf("%s\n",` > hypertorecdiff((1+x)^n*(1+y)^m/(1-x*y)^k,[n,m],f(x,y));`):

print((n+m+2-k)*f(n,m,x,y)+(-n-2+k)*f(n+1,m,x,y) = Diff(-x*(1+x)*f(n,m,x,y),x)+Diff((1+y)*f(n,m,x,y),y)):

print(` `):
printf("%s\n",`> hypertorecdiff((1+x)^n*(1+y)^m/(1-x*y)^k,[n,m],f(x,y),[0,1]);`):

print((n+m+2-k)*f(n,m,x,y)+(k-m-2)*f(n,m+1,x,y) = Diff((1+x)*f(n,m,x,y),x)+Diff(-y*(1+y)*f(n,m,x,y),y)):


print(` `):

print(` The following call asks MultInt to find a non-zero WZ-equation of with the given ansatz.`):

printf("%s\n",` > hypertorecdiff(1/(1-x-x^2-y-y^2)/x^(m+1)/y^(n+1),[n,m],f(x,y),[N,M],[N,M,N*M,M^2]);`): 

print((8+4*m+4*n)*f(n,m,x,y)+(2*n+2)*f(n+1,m,x,y)+(n+1)*f(n+1,m+1,x,y)+(8+2*n+4*m)*f(n,m+1,x,y)+(-10-5*m)*f(n,2+m,x,y) = Diff(-(4*x^2+4*x-5)*f(n,m,x,y)/x,x)+Diff(-(1+2*y)*(1+2*x)*f(n,m,x,y)/x,y)):

 
fi:


if nops([args])=1 and op(1,[args])=`checkrecdiff` then

print(` checkrecdiff: checks whether a  given WZ-equation satisfied by a given function. `):
print(` `):
printf("%s\n",`  Calling Sequence: `):
printf("%s\n",`    checkrecdiff(eq,f)`):
print(` `):
printf("%s\n",` Parameters:  `):
printf("%s\n",`    eq       - a recurrence-differential(WZ) equation`):
printf("%s\n",`    f        - algebraic expression   `):
print(` `):
printf("%s\n",`Description:  `):
printf("%s\n",` checkrecdiff  is useful,  in particular, to verify the output of the procedure hypertorecdiff.`):
printf("%s\n",`Whenever there  is  a  conflict  between  the  function  name  checkrecdiff and another name used in the same session, use the form MultInt['checkrecdiff']. `):
  
print(` `):
printf("%s\n",` Example:`):
printf("%s\n",` > with(MultInt): `):
printf("%s\n",`> rec:=(n+1)^2*(n-1)^2*f(n,x,y)-3*(2+3*n)*n^2*(1+3*n)*f(n+1,x,y) = Diff(-(1+x)*(x*y^2*n^3+x*y^2*n^2-x*y^2*n-x*y^2-6*n*y^2+10*n^3*y^2-72*n^3-60*n^2-4*n^2*y^2-12*n)/x/y^2*f(n,x,y),x)+Diff(n*(1+y)*(12*x*y^2*n^2+22*n*y^2+6*y^2+2*x*y^2+20*n^2*y^2+10*x*y^2*n-39*n*y-33*n^2*y-9*n^2*x*y-3*n*x*y-12*y+3*n+3*n*x+9*x*n^2+9*n^2)/x^2/y^2*f(n,x,y),y);`):

rec:=(n+1)^2*(n-1)^2*f(n,x,y)-3*(2+3*n)*n^2*(1+3*n)*f(n+1,x,y) = Diff(-(1+x)*(x*y^2*n^3+x*y^2*n^2-x*y^2*n-x*y^2-6*n*y^2+10*n^3*y^2-72*n^3-60*n^2-4*n^2*y^2-12*n)/x/y^2*f(n,x,y),x)+Diff(n*(1+y)*(12*x*y^2*n^2+22*n*y^2+6*y^2+2*x*y^2+20*n^2*y^2+10*x*y^2*n-39*n*y-33*n^2*y-9*n^2*x*y-3*n*x*y-12*y+3*n+3*n*x+9*x*n^2+9*n^2)/x^2/y^2*f(n,x,y),y);

printf("%s\n",`> checkrecdiff(rec, (1+1/x)^n*(1+1/y)^(2*n)/x^n/y^n); `):
print(true):

fi:

if nops([args])=1 and op(1,[args])=`sumtointn` then

print(` sumtointn: finds an integral representation of a given binomial sum`):

print(` `):
printf("%s\n",` Calling Sequences: `):
printf("%s\n",`    sumtointn(binocoeffs, x, k, [l,u]) `):
printf("%s\n",`    sumtointn(binocoeffs, x, [k1,k2,...,km], [[l1,u1],...,[lm,um]]) `):
 
print(` `):
printf("%s\n",`  Parameters: `):
printf("%s\n",`    binocoeffs    - finite product of binomial coeffs `):
printf("%s\n",`    x             - algebraic expression `):
printf("%s\n",`    k, k1,...,km  - names, the summation variables`):
printf("%s\n",`    l,u,l1,u1,... - expressions, lower and upper bounds of each of the summation variables`):

print(` `):
           
printf("%s\n",` Description:`):

print(` `):

printf("%s\n",` sumtointn finds an integral representation  of a summation expression of the form:`):
print(Sum(`(product of binocoeffs) x`,k)):
printf("%s\n",` i.e, finds a constant term (CT) expression: CT(F(z1,...,zr)) of the given sum and ouputs F with the integration variables z1,...,zr.`):
printf("%s\n",`Whenever there is a conflict between the function name  msumtointn  and another name used in the same session, use the long form MultInt['sumtointn']`):

printf("%s\n",` Examples: `):

printf("%s\n",`To find an intergral representation of:`):
print(Sum((-1)^k*binomial(n,k)^3, k)):

printf("%s\n",`> with(MultInt): `):
printf("%s\n",`> sumtointn(binomial(n,k)^3,(-1)^k,k,[0,infinity]);`):   
print(((z_1+1)/z_1)^n*((z_2+1)/z_2)^n*(1-z_1*z_2)^n, [z_1, z_2]):

printf("%s\n",`To find a  CT representaion of :`):

print(Sum((-1)^(i+j)*binomial(i+j,i)*binomial(n,i)*binomial(n,j),`i,j`)):

printf( "%s\n",`> with(MultInt):`):
printf("%s\n",`> sumtointn(binomial(i+j,i)*binomial(n,i)*binomial(n,j),(-1)^(i+j),[i,j],[[0,infinity],[0,infinity]]);`):
print((-z_1)^n*(-1/z_1)^n, [z_1]):

fi:


if nops([args])=1 and op(1,[args])=`sumtorecdiff` then

print(`sumtorecdiff - finds a non-zero recurrence-differential (WZ) equation for  the Constant Term expression of a given sum.`):
print(` `):
printf("%s\n", `Calling Sequences:  `): 
printf("%s\n", `  sumtorecdiff(binocoeffs, x, k, n, [l,u]) `):
printf("%s\n", `  sumtorecdiff(binocoeffs, x, [k1,...,km], n, [[l1,u1],...,[lm,um]])`):
printf("%s\n", `  sumtorecdiff(binocoeffs, x, k, n, [l,u], rec_order) `):
printf("%s\n", `  sumtorecdiff(binocoeffs, x, [k1,...,km], n, [[l1,u1],...,[lm,um]], rec_order)`):

print(` `):
printf("%s\n", ` Parameters:  `):
printf("%s\n", `    binocoeffs   - finite product of binomial coefficients`):
printf("%s\n", `    x            - algebraic expression `):
printf("%s\n", `    k,k1,...,km  - names, the summation variables`):
printf("%s\n", `    n            - name, the recurrence variable`):
printf("%s\n", `    l,u,l1,u1,... - expressions, lower and upper bounds of the summation variables`):
printf("%s\n", `    rec_order    - non-negative integer, the order of the recurrence`):

print(` `):
printf("%s\n", ` Description:`):
printf("%s\n", `    Given a  summation expression of the form `):

print(S(n) = Sum(Product(binocoeffs)*x,k)):

printf("%s\n", `sumtorecdiff  first finds an intergal representation  (constant term expression) of  the given sum : S(n) = CT(f(n,z1,...,zr)) and then finds a non-zero WZ-equation statisfied by f(n,z1,..,zr) and outputs the WZ-equation and f(n,z1,..,zn).`):

printf("%s\n",` - Whenever there  is  a  conflict  between  the  function  name  msum and another name used in the same session, use the long form  MultInt['sumtorecdiff']. `):

printf("%s\n", `Examples: `):
printf("%s\n",`To find a non-zero WZ-equation satisfied by the CT expresion of Dixon's sum`):
 print(Sum(binomial(2*n,k)^3*(-1)^k,k)):
printf("%s\n",`> with(MultInt):`):
printf("%s\n",`> sumtorecdiff(binomial(2*n,k)^3,(-1)^k,k,n,[0,infinity],1);`):

print([3*(3*n+1)*(2+3*n)*_F(n,z_1,z_2)+(n+1)^2*_F(n+1,z_1,z_2) = Diff(_R(z_1,z_2)*_F(n,z_1,z_2),z_1)+Diff(_R(z_2,z_1)*_F(n,z_1,z_2),z_2), _R(z_1,z_2) = -1/4*(z_1+1)*(1+3*n-4*z_1^3*z_2^3+4*z_1^2*z_2^3*n+4*z_1^2*z_2^3*n^2+42*z_1^2*z_2^2*n+7*z_1+9*z_1^3*z_2*n-4*z_1^3*z_2*n^2-12*z_1^3*z_2*n^3+4*z_1^2*z_2*n^2-9*z_1^2*n*z_2+12*z_1^2*z_2*n^3-10*z_1*z_2^2-34*z_1*n*z_2^2-4*z_1^2*z_2+4*z_1^3*z_2+4*z_1*z_2+10*z_1^2*z_2^2-z_1^3*z_2^4+2*n^2+8*z_1*z_2^4+21*z_1*n+12*z_1*z_2^2*n^3-32*z_1*z_2^4*n^3+10*z_1*z_2^4*n^2-24*z_1*z_2^4*n^4+26*z_1*z_2^4*n-12*n^3*z_1*z_2-4*n^2*z_1*z_2-2*z_1^3*z_2^4*n^2-18*z_1*z_2^2*n^2+12*z_2^3*n^3+2*z_2^2-4*z_2^3+z_1^2*z_2^4+24*z_1^2*z_2^4*n^2+12*z_1^2*z_2^4*n^3+13*z_1^2*z_2^4*n-18*z_1^3*z_2^3*n-8*z_1^2-2*z_1^3*z_2^2-3*n*z_2+2*z_1*n^2+22*z_2^2*n^2-3*z_1^3*z_2^4*n-12*n^3*z_1^2*z_2^2-24*n^4*z_1^2*z_2^2-30*z_1*z_2^3*n^2+12*z_1*z_2^3*n^3+24*z_1*z_2^3*n^4-25*z_1*z_2^3*n-22*z_1^3*z_2^2*n^2+8*z_1^3*z_2^2*n^3+24*z_1^3*z_2^2*n^4-14*z_1^3*z_2^2*n-26*z_1^3*z_2^3*n^2-12*z_1^3*z_2^3*n^3-24*z_1*n^4+24*n^3*z_2+44*n^2*z_1^2*z_2^2-4*z_1*z_2^3-36*z_1*n^3+9*n*z_1*z_2-10*z_1^2*n^2+24*n^4*z_2+14*z_2^2*n-8*z_2^2*n^3-24*z_2^2*n^4+4*z_2^3*n^2-9*z_2^3*n+32*z_1^2*n^3+24*z_1^2*n^4-26*z_1^2*n)/(z_2^2*(2*n+1)*z_1), _F(n,z_1,z_2) = -1/4*((z_1+1)/z_1)^(2*n)*((z_2+1)/z_2)^(2*n)*(1-z_1*z_2)^(2*n)/(z_1*z_2*Pi^2)]):

printf( "%s\n",`To find a non-zero WZ-equation for the CT form of :`):
print(Sum(binomial(i+j,i)*binomial(n,i)*binomial(n,j)*(-1)^(i+j),`i,j`)):

printf( "%s\n",`> with(MultInt):`): 
printf( "%s\n",`> sumtorecdiff(binomial(i+j,i)*binomial(n,i)*binomial(n,j),(-1)^(i+j),[i,j],n,[[0,infinity],[0,infinity]]);`):
print([_F(n,z_1)-_F(n+1,z_1) = Diff(0,z_1), _F(n,z_1) = -1/2*I*(-z_1)^n*(-1/z_1)^n/(z_1*Pi)]):
fi:


if nops([args])=1 and op(1,[args])=`esp` then

print(`  esp: generates the elementary symmetric polynomials of a given order`):
print( `  `):
printf( "%s\n",` Calling Sequence: esp(n,[v1,...,vr])`):

print(` `):

printf( "%s\n",` Parameters:  `):
printf( "%s\n",`    n              - a positive integer`):
printf( "%s\n",`    v1,v2,...,vr   - names  `):
print(` `):
printf ("%s\n",` Description: `):
printf( "%s\n",`- The function esp(n,[v1,...,vr]) outputs the elementary symmetric polynomial in the variables [v1,...,vr] of order n.`):
printf( "%s\n",`- Whenever there is a  conflict  between  the  function  name  esp  and another name used in the same session, use the long form MultInt['esp']. `):
 
print(` `):
printf( "%s\n",`Examples: `):
printf( "%s\n",`> with(MultInt):`):
printf( "%s\n",`> esp(2,[x,y,z]);`):
print(x*y+x*z+y*z):
fi:


if nops([args])=1 and op(1,[args])=`sym` then

print(` sym: checks whether a given function is symmetric w.r.t. given variables.`):
print(` `):
printf("%s\n",` Calling  Sequence: `):
printf("%s\n",`     sym(f,[v1,v2,...,vn])`):
print(` `):
printf("%s\n",` Parameters:  `):

printf("%s\n",`    f           - an algebraic expression`):
prinf("%s\n",`    v1,...,vn    - names`):

print(` `):
   
printf("%s\n",` Description: `):

printf("%s\n",` - sym returns true,  if  f  is symmetric w.r.t v1,..,vn, false otherwise. `):
printf("%s\n",` - Whenever there  is  a  conflict  between  the  function  name  sym  and another name used in the same session, use the form MultInt['sym']. `):

print(` `):

printf("%s\n",` Example: `):
printf("%s\n",`> with(MultInt):`):
printf("%s\n",`> sym((1-1/x)^n*(1-1/y)^n/x/y,[x,y]) ;`):
print(true):
fi:

end: