#file MultInt.mpl #last update: Dec, 2000. # 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 # Temple University, Math Dept. # akalu@math.temple.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 6. # maple.hdb is the accompanying HELP file. After downloading it, it will be accessible # via Maple by adding the directory path where this file is stored to the libname. # Example: In a Maple session type # # libname; # # to see its current value. # # # Suppose you stored the file maple.hdb in the directory /home/subdir1/. Then # add this path to the libname as follows: # # libname:= libname, `/home/subdir1`; # # After executing this command to see, for instance, the help page of MultInt type: # # ?MultInt # # in your Maple session. # 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=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 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(_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(`solve/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(`solve/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(`solve/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(`solve/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: 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(``): 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(``): RETURN(ope=sum('(Diff(op(i,Rlst)*fnam(mainvar,op(intnvars)),intnvars[i]))', 'i'=1..nops(intnvars))): 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: 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(``): 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.`): 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: 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(``): 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(``): 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)])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:`): print(` ?MultInt or help(MultInt) `): print(` `): print(` For help with a specific procedure, type:`): print(` ?procedure name or help(procedure name) `): print(``): print(` CAUTION: this version of MultInt is for Maple 6.`): print(``): if not member(`.`,[libname]) then libname:=libname,`.`: fi: