/* FPS 1.0 - A MuPAD-Package for Computations of Formal Power Series (C) Torsten Sprenger, University of Kassel, 2006 Torsten Sprenger sprenger@mathematik.uni-kassel.de http://www.mathematik.uni-kassel.de/~sprenger */ /*************************************************************************************************/ pochhammer:=proc(n,k)/* defines the pochhammer symbol */ local j; begin if testtype(k,Type::NonNegInt) then return(_mult(n+j $ j=0..k-1)) elif testtype(n,Type::PosInt) then return((k+n-1)!/(n-1)!) elif testtype(n,Type::Rational) then if denom(n)=2 and testtype(n,Type::Positive) then if n=1/2 then return((2*k)!/(4^k*k!)) else return(simplify(expand(pochhammer(1/2,_div(numer(n),2)+k)/pochhammer(1/2,_div(numer(n),2))))) end_if; else return(procname(args())) end_if; else return(procname(args())) end_if; end_proc: DEtoPol:=proc(De,Fx,y)/* converts a differential equation in a polynomial */ local F,x,i,n; begin F:=op(Fx,0); x:=op(Fx,1); n:=ode::getOrder(De,F(x)); subsex(De,diff(F(x),x$n-i)=y^(n-i) $ i=0..n) end_proc: PoltoDE:=proc(Pol,y,Fx)/* converts a polynomial in a differential equation */ local F,x,i,n; begin F:=op(Fx,0); x:=op(Fx,1); n:=degree(Pol,[y]); subsex(Pol,y^(n-i)=diff(F(x),x$n-i) $ i=0..n) end_proc: ReducePol:=proc(Pol,x,redPol)/* reduces Pol in x by the rule of redPol */ local redPOL,redlist,POL,n,r; begin redPOL:=poly(redPol/lcoeff(redPol,[x]),[x]); n:=degree(redPOL); POL:=poly(Pol,[x]); redlist:=map(poly2list(-redPOL+lmonomial(redPOL)),m->[m[1],n-m[2]]); while degree(POL)>=n do list:=poly2list(POL); list:=map(list,m->_if(m[2]>=n,([m[1]*r[1], m[2]-r[2]] $ r in redlist),m)); POL:=poly(list,[x]); end_while; expr(POL); end_proc: ReducePol:=proc(Pol,x,redPol)/* reduces Pol in x by the rule of redPol */ local m,n,redPOL; begin n:=degree(redPol,[x]); m:=degree(Pol,[x]); redPOL:=-(redPol-lmonomial(redPol,[x])); if m>=n then ReducePol(subs(expand(Pol),x^m=x^(m-n)*redPOL),x,redPol) else collect(Pol,x,factor) end_if; end_proc: ReduceDE:=proc(De,Fx,redDe)/* reduces De in Fx by the rule of redDE */ local F,x,m,n,i,temp,redDE,DE; begin F:=op(Fx,0); x:=op(Fx,1); DE:=expand(De); n:=ode::getOrder(redDe,F(x)); m:=ode::getOrder(DE,F(x)); temp:=genident(); redDE:=op(solve(subsex(redDe,diff(F(x),x$n)=temp),temp,IgnoreSpecialCases)); if m>=n then return(ReduceDE(subs(DE,diff(F(x),x$m)=diff(redDE,x$m-n)),F(x),redDe)) else return(collect(DE,[diff(F(x),x$i) $ i=0..n-1],factor)) end_if; end_proc: ReduceRE:=proc(Re,Ak,redRe)/* reduces Re in Ak by the rule of redRE */ local A,k,m,n,i,redRE,RE; begin A:=op(Ak,0); k:=op(Ak,1); RE:=expand(Re); n:=OrderRE(redRe,A(k)); m:=OrderRE(RE,A(k)); redRE:=op(solve(redRe,A(k+n),IgnoreSpecialCases)); if m>=n then return(ReduceRE(subsex(RE,A(k+m)=subs(redRE,k=k+m-n)),A(k),redRe)) else return(collect(RE,[A(k+i) $ i=0..n-1],factor)) end_if; end_proc: FallFact:=proc(a,k)/* determines coefficients (w.r.t. falling factorials) of a polynomial a in k */ local pola,n,coe,coelist; begin pola:=expand(a); n:=degree(pola,k); coelist:=[]; for i from 0 to n do coe:=coeff(pola,k,n-i); pola:=expand(pola-coe*pochhammer(k-(n-i)+1,n-i)); coelist:=[op(coelist),coe]; end_for; coelist end_proc: SUM:=newDomain("SUM"):/* output of FPS */ SUM::new:=(s1,s2)->new(SUM,s1,s2): SUM::print:=(s)->hold(sum)(op(s,1),op(s,2)): SUM::term:=(s)->op(s,1): /*************************************************************************************************/ /* Global Variables */ MAXDEORDER:=5: MAXREORDER:=5: /*************************************************************************************************/ /* HolonomicDE determines a differential equation for a holonomic function f */ HolonomicDE:=proc(f,Fx) local F,x,rat,a,n,i,j,k,DElist,filled,d,eq,sol,DE; begin F:=op(Fx,0); x:=op(Fx,1); // DEorder=1 rat:=normal(diff(f,x)/f); if testtype(rat,Type::RatExpr(x)) then return(factor(denom(rat))*diff(F(x),x)-factor(numer(rat))*F(x)=0) end_if; // DEorder>1 a:=[genident() $ i=1..MAXDEORDER+1]; if testtype(f,"_plus") then DElist:=[{a[1]*op(f,k)} $ k=1..nops(f)] // initialize DElist with f else DElist:=[{a[1]*f}] end_if; n:=1; while n<=MAXDEORDER do d:=diff(f,x$n); if not(testtype(d,"_plus")) then d:={d} end_if; for k from 1 to nops(d) do // checking l.d. for every summand of the n-th derivative filled:=FALSE; for i from 1 to nops(DElist) do if testtype(normal(op(d,k)/DElist[i][1]),Type::RatExpr(x)) then DElist[i]:={op(DElist[i]),a[n+1]*op(d,k)}; filled:=TRUE; end_if; end_for; if not(filled) then DElist:=[op(DElist),{a[n+1]*op(d,k)}] // construct a new equivalence class if no summand of the n-th derivative is l.d. end_if; end_for; eq:=subs(map(DElist,y->_plus(op(y))),a[n+1]=1); // add terms in equivalence classes sol:=solve(eq,{a[j] $ j=1..n},IgnoreSpecialCases); if nops(sol)<>0 then DE:=subs(_plus(a[j]*diff(F(x),x$j-1) $ j=1..n+1),op(sol)); DE:=collect(numer(normal(DE)),[diff(F(x),x$j-1) $ j=1..n+1],factor); DE:=subs(DE,{a[j]=1 $ j=1..n+1}); return(DE=0); end_if; n:=n+1; end_while; error("There exists no differential equation with order less or equal than MAXDEORDER= ".MAXDEORDER); end_proc: /*************************************************************************************************/ /* HolonomicDE determines a differential equation for a holonomic function f based on the product and sum algorithms ProductDE and SumDE */ HolonomicDE_:=proc(f,Fx) local F,x,rat,a,n,i,j,k,DElist,filled,d,eq,sol,DE,baseDE; begin F:=op(Fx,0); x:=op(Fx,1); // Product if testtype(f,"_mult") then DE:=HolonomicDE_(op(f,1),F(x)); for i from 2 to nops(f) do DE:=ProductDE(DE,HolonomicDE_(op(f,i),F(x)),F(x)) end_for; return(DE); end_if; // Sum if testtype(f,"_plus") then DE:=HolonomicDE_(op(f,1),F(x)); for i from 2 to nops(f) do DE:=SumDE(DE,HolonomicDE_(op(f,i),F(x)),F(x)) end_for; return(DE); end_if; // Power if testtype(f,"_power") and testtype(op(f,2),Type::PosInt) then baseDE:=HolonomicDE_(op(f,1),F(x)); DE:=baseDE; for i from 2 to op(f,2) do DE:=ProductDE(DE,baseDE,F(x)) end_for; return(DE); end_if; // DEorder=1 rat:=normal(diff(f,x)/f); if testtype(rat,Type::RatExpr(x)) then return(factor(denom(rat))*diff(F(x),x)-factor(numer(rat))*F(x)=0) end_if; // DEorder>1 a:=[genident() $ i=1..MAXDEORDER+1]; DElist:=[{a[1]*f}]; // initialize DElist with f n:=1; while n<=MAXDEORDER do d:=diff(f,x$n); if not(testtype(d,"_plus")) then d:={d} end_if; for k from 1 to nops(d) do // checking l.d. for every summand of the n-th derivative filled:=FALSE; for i from 1 to nops(DElist) do if testtype(normal(op(d,k)/DElist[i][1]),Type::RatExpr(x)) then DElist[i]:={op(DElist[i]),a[n+1]*op(d,k)}; filled:=TRUE; end_if; end_for; if not(filled) then DElist:=[op(DElist),{a[n+1]*op(d,k)}] // construct a new equivalence class if no summand of the n-th derivative is l.d. end_if; end_for; eq:=subs(map(DElist,y->_plus(op(y))),a[n+1]=1); // add terms in equivalence classes sol:=solve(eq,{a[j] $ j=1..n},IgnoreSpecialCases); if nops(sol)<>0 then DE:=subs(_plus(a[j]*diff(F(x),x$j-1) $ j=1..n+1),op(sol)); DE:=collect(numer(normal(DE)),[diff(F(x),x$j-1) $ j=1..n+1],factor); DE:=subs(DE,{a[j]=1 $ j=1..n+1}); return(DE=0); end_if; n:=n+1; end_while; error("There exists no differential equation with order less or equal than MAXDEORDER= ".MAXDEORDER); end_proc: /*************************************************************************************************/ /* DEtoRE converts a differential equation De of a formal power series Fx in a recurrence equation for the coefficients Ak */ DEtoRE:=proc(De,Fx,Ak) local F,x,A,k,RE,DE,i,j,m,mn,mx; begin F:=op(Fx,0); x:=op(Fx,1); A:=op(Ak,0); k:=op(Ak,1); DE:=expand(lhs(De)-rhs(De)); RE:=0; mn:=infinity; mx:=-infinity; for i from 1 to nops(DE) do j:=degree(op(DE,i),x); m:=ode::getOrder(op(DE,i),F(x)); RE:=RE+subsex(op(DE,i),x^j*diff(F(x),x$m)=pochhammer(k+1-j,m)*A(k+m-j)); // apply rule: x^j F^(m) -> pochhammer(k+1-j,m)*A(k+m-j) mn:=min(mn,m-j); mx:=max(mx,m-j); end_for; RE:=subs(RE,k=k-mn); collect(RE,[A(k+j) $ j=0..mx-mn],factor)=0 end_proc: /*************************************************************************************************/ /* REtoDE converts a recurrence equation Re for the coefficients Ak in a differential equation for the formal power series Fx */ REtoDE:=proc(Re,Ak,Fx) local F,x,A,k,RE,DE,coe,i,j,m,l; begin F:=op(Fx,0); x:=op(Fx,1); A:=op(Ak,0); k:=op(Ak,1); RE:=expand(lhs(Re)-rhs(Re)); DE:=0; for i from 1 to nops(RE) do j:=degree(op(RE,i),k); m:=OrderRE(op(RE,i),A(k)); coe:=FallFact((k-m)^j,k); // (k-m)^j = sum(b_l fallfact(k,l),l=0..j) DE:=DE+subsex(op(RE,i),k^j*A(k+m)=_plus(coe[j-l+1]*x^(l-m)*diff(F(x),x$l) $ l=0..j)); // apply rule: k^j A(k+m) -> sum(b_l x^(l-m) F^(l),l=0..j) end_for; DE:=numer(normal(DE)); collect(DE,[diff(F(x),x$j) $ j=0..degree(RE,k)],factor)=0 end_proc: /*************************************************************************************************/ /* SumtoHyper converts a hypergeometric series with summand a into the notation of a generalized hypergeometric function pFq */ SumtoHyper:=proc(a,k) local rat,num,den,low,upp,x,i,j,shift,minlow,minupp,coe1,coe0; begin rat:=simplify(expand(subs(a,k=k+1)/a)); num:=factor(numer(rat)); den:=factor(denom(rat)); x:=1; // determine upper parameters upp:=[];minupp:=infinity; if nops(num)>2 then x:=x*op(num,1) else x:=x*num; end_if; for i from 2 to nops(num) step 2 do coe1:=coeff(op(num,i),k,1); coe0:=coeff(op(num,i),k,0); if coe1<>0 then x:=x*coe1^op(num,i+1); if testtype(coe0/coe1,Type::Integer) then minupp:=min(minupp,coe0/coe1) end_if; for j from 1 to op(num,i+1) do upp:=[op(upp),coe0/coe1] end_for; else x:=x*op(num,i)^op(num,i+1) end_if; end_for; // determine lower parameters low:=[];minlow:=infinity; if nops(den)>2 then x:=x/op(den,1) else x:=x/den; end_if; for i from 2 to nops(den) step 2 do coe1:=coeff(op(den,i),k,1); coe0:=coeff(op(den,i),k,0); if coe1<>0 then x:=x/coe1^op(den,i+1); if testtype(coe0/coe1,Type::Integer) then minlow:=min(minlow,coe0/coe1) end_if; for j from 1 to op(den,i+1) do low:=[op(low),coe0/coe1] end_for; else x:=x/op(den,i)^op(den,i+1) end_if; end_for; // determine shift shift:=min(minlow,minupp); if testtype(shift,Type::NegInt) or shift=0 then shift:=shift-1; upp:=map(upp,y->y-shift); low:=map(low,y->y-shift); else shift:=0 end_if; // include or exclude 1 if has(low,1) then delete low[contains(low,1)] else upp:=[op(upp),1] end_if; Simplify(subs(a,k=-shift))*Hypergeom(upp,low,x) end_proc: /*************************************************************************************************/ /* SumDE determines a differential equation out of two differential equations De1 and De2, which is satisfied by the sum of the solutions of De1 and De2 */ /*SumDE:=proc(De1,De2,Fx) local F,x,G,H,i,l,m,n,eq,J,sol,DE,DE1,DE2,rules,sumDE,vars,varG,varH; begin F:=op(Fx,0); x:=op(Fx,1); m:=ode::getOrder(De1,F(x)); n:=ode::getOrder(De2,F(x)); G:=genident();H:=genident(); sumDE:=[genident() $ i=1..n+m+1]; varG:=[genident() $ i=1..m]; varH:=[genident() $ i=1..n]; rules:=diff(G(x),x$m-i)=varG[i] $ i=1..m,diff(H(x),x$n-i)=varH[i] $ i=1..n; vars:={op(varG),op(varH)}; J:=max(m,n); DE1:=subs(lhs(De1)-rhs(De1),F(x)=G(x)); DE2:=subs(lhs(De2)-rhs(De2),F(x)=H(x)); eq:={}; for i from 1 to J-1 do eq:={op(eq),subs(numer(normal(sumDE[i]- ReduceDE(diff(G(x),x$i-1),G(x),DE1)- ReduceDE(diff(H(x),x$i-1),H(x),DE2))),rules)} end_for; for i from J to n+m+1 do eq:={op(eq),subs(numer(normal(sumDE[i]- ReduceDE(diff(G(x),x$i-1),G(x),DE1)- ReduceDE(diff(H(x),x$i-1),H(x),DE2))),rules)}; sol:=groebner::eliminate(eq,vars); if nops(sol)<>0 then DE:=collect(op(sol),sumDE,factor); DE:=subs(DE,sumDE[l]=diff(F(x),x$l-1) $ l=1..i)=0; return(DE); end_if; end_for; end_proc:*/ SumDE:=proc(De1,De2,Fx) local F,x,G,H,i,j,m,n,eq,J,sol,DE,DE1,DE2,coe,vars,approach,LCM,dencoe,rat; begin F:=op(Fx,0); x:=op(Fx,1); m:=ode::getOrder(De1,F(x)); n:=ode::getOrder(De2,F(x)); G:=genident();H:=genident(); coe:=[genident() $ i=1..n+m+1]; vars:=[diff(G(x),x$i-1) $ i=1..m, diff(H(x),x$i-1) $ i=1..n]; J:=max(m,n); DE1:=subs(lhs(De1)-rhs(De1),F(x)=G(x)); DE2:=subs(lhs(De2)-rhs(De2),F(x)=H(x)); approach:=0; for i from 1 to J-1 do rat:=ReduceDE(diff(G(x),x$i-1),G(x),DE1)+ ReduceDE(diff(H(x),x$i-1),H(x),DE2); dencoe[i]:=denom(rat); approach:=approach+coe[i]*numer(rat); end_for; for i from J to n+m+1 do rat:=ReduceDE(diff(G(x),x$i-1),G(x),DE1)+ ReduceDE(diff(H(x),x$i-1),H(x),DE2); dencoe[i]:=denom(rat); approach:=approach+coe[i]*numer(rat); eq:=numer(normal(approach)); eq:=collect(eq,vars,factor); eq:={coeff(eq,vars)}; sol:=solve(eq,{coe[j] $ j=1..i},IgnoreSpecialCases); if subs({coe[j] $ j=1..i},op(sol))={0} then else LCM:=lcm(dencoe[j] $ j=1..i); DE:=_plus(coe[j]*dencoe[j]/LCM*diff(F(x),x$j-1) $ j=1..i); DE:=subs(DE,op(sol)); DE:=subs(DE,map(indets(DE) minus _union(indets(De1),indets(De2)),z->z=1)); DE:=numer(normal(DE)); DE:=collect(DE,[diff(F(x),x$j-1) $ j=1..i],factor)=0; return(DE); end_if; end_for; end_proc: /*************************************************************************************************/ /* ProductDE determines a differential equation out of two differential equations De1 and De2, which is satisfied by the product of the solutions of De1 and De2 */ /*ProductDE:=proc(De1,De2,Fx) local F,x,G,H,i,j,l,m,n,eq,J,sol,DE,DE1,DE2,rules,productDE,vars,redDE; begin F:=op(Fx,0); x:=op(Fx,1); m:=ode::getOrder(De1,F(x)); n:=ode::getOrder(De2,F(x)); G:=genident();H:=genident(); productDE:=[genident() $ i=1..n*m+1]; vars:=[genident() $ i=1..n*m]; rules:=diff(G(x),x$m-i)*diff(H(x),x$n-j)=vars[m*(j-1)+i] $ i=1..m $ j=1..n; J:=max(m,n); DE1:=subs(lhs(De1)-rhs(De1),F(x)=G(x)); DE2:=subs(lhs(De2)-rhs(De2),F(x)=H(x)); eq:={}; redDE:=G(x)*H(x); for i from 1 to J-1 do eq:={op(eq),subsex(numer(normal(productDE[i]-expand(redDE))),rules)}; redDE:=diff(redDE,x); redDE:=ReduceDE(redDE,G(x),DE1); redDE:=ReduceDE(redDE,H(x),DE2); end_for; for i from J to n*m+1 do eq:={op(eq),subsex(numer(normal(productDE[i]-expand(redDE))),rules)}; sol:=groebner::eliminate(eq,vars); if nops(sol)<>0 then DE:=collect(op(sol),productDE,factor); DE:=subs(DE,productDE[l]=diff(F(x),x$l-1) $ l=1..i)=0; return(DE); end_if; redDE:=diff(redDE,x); redDE:=ReduceDE(redDE,G(x),DE1); redDE:=ReduceDE(redDE,H(x),DE2); end_for; end_proc:*/ ProductDE:=proc(De1,De2,Fx) local F,x,G,H,i,j,m,n,eq,J,sol,DE,DE1,DE2,coe,dencoe,vars,redDE,rules,approach,LCM,rat; begin F:=op(Fx,0); x:=op(Fx,1); m:=ode::getOrder(De1,F(x)); n:=ode::getOrder(De2,F(x)); G:=genident();H:=genident(); coe:=[genident() $ i=1..n*m+1]; vars:=[genident() $ i=1..n*m]; rules:=diff(G(x),x$m-i)*diff(H(x),x$n-j)=vars[m*(j-1)+i] $ i=1..m $ j=1..n; J:=max(m,n); DE1:=subs(lhs(De1)-rhs(De1),F(x)=G(x)); DE2:=subs(lhs(De2)-rhs(De2),F(x)=H(x)); approach:=0; redDE:=G(x)*H(x); for i from 1 to J-1 do rat:=normal(expand(redDE)); dencoe[i]:=denom(rat); approach:=approach+subsex(coe[i]*numer(rat),rules); redDE:=diff(redDE,x); redDE:=ReduceDE(redDE,G(x),DE1); redDE:=ReduceDE(redDE,H(x),DE2); end_for; for i from J to n*m+1 do rat:=normal(expand(redDE)); dencoe[i]:=denom(rat); approach:=approach+subsex(coe[i]*numer(rat),rules); eq:=numer(normal(approach)); eq:=collect(eq,vars,factor); eq:={coeff(eq,vars)}; sol:=solve(eq,{coe[j] $ j=1..i},IgnoreSpecialCases); if subs({coe[j] $ j=1..i},op(sol))={0} then else LCM:=lcm(dencoe[j] $ j=1..i); DE:=_plus(coe[j]*dencoe[j]/LCM*diff(F(x),x$j-1) $ j=1..i); DE:=subs(DE,op(sol)); DE:=subs(DE,map(indets(DE) minus _union(indets(De1),indets(De2)),z->z=1)); DE:=numer(normal(DE)); DE:=collect(DE,[diff(F(x),x$j-1) $ j=1..i],factor)=0; return(DE); end_if; redDE:=diff(redDE,x); redDE:=ReduceDE(redDE,G(x),DE1); redDE:=ReduceDE(redDE,H(x),DE2); end_for; end_proc: /*************************************************************************************************/ /* AlgebraicDE determines a differential equation for an algebraic function Fx out of an equation eq containing Fx and x */ AlgebraicDE:=proc(eq,Fx) local F,x,expr,i,j,k,n,temp,der,a,rules,approach,pol,sol,DE; begin F:=op(Fx,0); x:=op(Fx,1); if has(eq,F(x)) then expr:=eq else expr:=subs(eq,F=F(x)) end_if; n:=degree(expr,[F(x)]); a:=[genident() $ i=1..n+1]; rules:=[]; temp:=genident(); for i from 1 to n do der:=op(solve(subsex(diff(expr,x$i),diff(F(x),x$i)=temp),temp,IgnoreSpecialCases)); der:=normal(subs(der,op(rules))); rules:=[diff(F(x),x$i)=der,op(rules)]; approach:=_plus(a[k]*diff(F(x),x$k-1) $k=1..i)+diff(F(x),x$i); pol:=numer(normal(subs(approach,op(rules)))); pol:=expand(divide(pol,expand(expr),[F(x)])[2]); // reduces pol sol:=solve({coeff(pol,[F(x)])},{a[j] $ j=1..i+1},IgnoreSpecialCases); if nops(sol)<>0 then DE:=subs(_plus(a[j]*diff(F(x),x$j-1) $ j=1..i+1),op(sol)); DE:=collect(numer(normal(DE)),[diff(F(x),x$j-1) $ j=1..i+1],factor); DE:=subs(DE,{a[j]=1 $ j=1..i+1}); return(DE=0); end_if; end_for; end_proc: /*************************************************************************************************/ /* SubsDE substitutes function f in a differential equation De for Fx */ SubsDE:=proc(De,Fx,f) local F,x,i,j,DE,n,rules,sol; begin F:=op(Fx,0); x:=op(Fx,1); DE:=lhs(De)-rhs(De); rules:={F(f)=F(x)}; n:=ode::getOrder(DE,F(x)); DE:=subs(DE,x=f); for i from 1 to n do sol:=solve(diff(F(f),x$i)=diff(F(x),x$i),(D@@i)(F)(f),IgnoreSpecialCases); rules:=map({(D@@i)(F)(f)=op(sol),op(rules)},y->lhs(y)=normal(subs(rhs(y),rules))); end_for; rules:=subs(rules,{(D@@i)(F)(f)=diff(F(f),f$i) $ i=1..n}); DE:=numer(normal(subs(DE,rules))); DE:=collect(DE,[diff(F(x),x$j) $ j=0..n],factor)=0; end_proc: /*************************************************************************************************/ /* HolonomicRE determines a recurrence equation for a holonomic series f */ HolonomicRE:=proc(f,Ak) local A,k,rat,a,b,n,i,j,RElist,filled,r,eq,sol,RE; begin A:=op(Ak,0); k:=op(Ak,1); a:=expand(f); // REorder=1 rat:=simplify(expand(subs(a,k=k+1)/a)); if testtype(rat,Type::RatExpr(k)) then return(factor(denom(rat))*A(k+1)-factor(numer(rat))*A(k)=0) end_if; // REorder>1 b:=[genident() $ i=1..MAXREORDER+1]; if testtype(a,"_plus") then RElist:=[{b[1]*op(a,j)} $ j=1..nops(a)] else RElist:=[{b[1]*a}] end_if; n:=1; while n<=MAXREORDER do r:=subs(a,k=k+n); if not(testtype(r,"_plus")) then r:={r} end_if; for j from 1 to nops(r) do filled:=FALSE; for i from 1 to nops(RElist) do if testtype(simplify(expand(op(r,j)/RElist[i][1])),Type::RatExpr(k)) then RElist[i]:={op(RElist[i]),b[n+1]*op(r,j)}; filled:=TRUE; end_if; end_for; if not(filled) then RElist:=[op(RElist),{b[n+1]*op(r,j)}] end_if; end_for; eq:=subs(map(RElist,y->_plus(op(y))),b[n+1]=1); sol:=solve(eq,{b[j] $ j=1..n},IgnoreSpecialCases); if nops(sol)<>0 then RE:=subs(_plus(b[j]*A(k+j-1) $ j=1..n+1),op(sol)); RE:=collect(RE,[A(k+j-1) $ j=1..n+1],x->simplify(expand(x))); RE:=collect(numer(normal(RE)),[A(k+j-1) $ j=1..n+1],factor); RE:=subs(RE,{b[j]=1 $ j=1..n+1}); return(RE=0); end_if; n:=n+1; end_while; error("There exists no recurrence equation with order less or equal than MAXREORDER= ".MAXREORDER); end_proc: /*************************************************************************************************/ /* HolonomicRE determines a recurrence equation for a holonomic function f based on the product and sum algorithms ProductRE and SumRE */ HolonomicRE_:=proc(f,Ak) local A,k,a,rat,b,n,i,j,RElist,filled,r,eq,sol,RE,baseRE; begin A:=op(Ak,0); k:=op(Ak,1); a:=f; // Product if testtype(a,"_mult") then RE:=HolonomicRE_(op(a,1),A(k)); for i from 2 to nops(a) do RE:=ProductRE(RE,HolonomicRE_(op(a,i),A(k)),A(k)) end_for; return(RE); end_if; // Sum if testtype(a,"_plus") then RE:=HolonomicRE_(op(a,1),A(k)); for i from 2 to nops(a) do RE:=SumRE(RE,HolonomicRE_(op(a,i),A(k)),A(k)) end_for; return(RE); end_if; // Power if testtype(a,"_power") and testtype(op(a,2),Type::PosInt) then baseRE:=HolonomicRE_(op(a,1),A(k)); RE:=baseRE; for i from 2 to op(a,2) do RE:=ProductRE(RE,baseRE,A(k)) end_for; return(RE); end_if; // REorder=1 rat:=simplify(expand(subs(a,k=k+1)/a)); if testtype(rat,Type::RatExpr(k)) then return(factor(denom(rat))*A(k+1)-factor(numer(rat))*A(k)=0) end_if; // REorder>1 b:=[genident() $ i=1..MAXREORDER+1]; RElist:=[{b[1]*a}]; n:=1; while n<=MAXREORDER do r:=subs(a,k=k+n); if not(testtype(r,"_plus")) then r:={r} end_if; for j from 1 to nops(r) do filled:=FALSE; for i from 1 to nops(RElist) do if testtype(simplify(expand(op(r,j)/RElist[i][1])),Type::RatExpr(k)) then RElist[i]:={op(RElist[i]),b[n+1]*op(r,j)}; filled:=TRUE; end_if; end_for; if not(filled) then RElist:=[op(RElist),{b[n+1]*op(r,j)}] end_if; end_for; eq:=subs(map(RElist,y->_plus(op(y))),b[n+1]=1); sol:=solve(eq,{b[j] $ j=1..n},IgnoreSpecialCases); if nops(sol)<>0 then RE:=subs(_plus(b[j]*A(k+j-1) $ j=1..n+1),op(sol)); RE:=collect(numer(normal(RE)),[A(k+j-1) $ j=1..n+1],factor); RE:=subs(RE,{b[j]=1 $ j=1..n+1}); return(RE=0); end_if; n:=n+1; end_while; error("There exists no recurrence equation with order less or equal than MAXREORDER= ".MAXREORDER); end_proc: /*************************************************************************************************/ /* SumRE determines a recurrence equation out of two recurrence equations Re1 and Re2, which is satisfied by the sum of the solutions of Re1 and Re2 */ /*SumRE:=proc(Re1,Re2,Ak) local A,k,B,C,i,l,m,n,eq,J,sol,RE,RE1,RE2,rules,sumRE,vars,varB,varC; begin A:=op(Ak,0); k:=op(Ak,1); m:=OrderRE(Re1,A(k)); n:=OrderRE(Re2,A(k)); B:=genident();C:=genident(); sumRE:=[genident() $ i=1..n+m+1]; varB:=[genident() $ i=1..m]; varC:=[genident() $ i=1..n]; rules:=B(k+i-1)=varB[i] $ i=1..m,C(k+i-1)=varC[i] $ i=1..n; vars:={op(varB),op(varC)}; J:=max(m,n); RE1:=subs(lhs(Re1)-rhs(Re1),A=B); RE2:=subs(lhs(Re2)-rhs(Re2),A=C); eq:={}; for i from 1 to J-1 do eq:={op(eq),subs(numer(normal(sumRE[i]- ReduceRE(B(k+i-1),B(k),RE1)- ReduceRE(C(k+i-1),C(k),RE2))),rules)} end_for; for i from J to n+m+1 do eq:={op(eq),subs(numer(normal(sumRE[i]- ReduceRE(B(k+i-1),B(k),RE1)- ReduceRE(C(k+i-1),C(k),RE2))),rules)}; sol:=groebner::eliminate(eq,vars); if nops(sol)<>0 then RE:=collect(op(sol),sumRE,factor); RE:=subs(RE,sumRE[l]=A(k+l-1) $ l=1..i)=0; return(RE); end_if; end_for; end_proc:*/ SumRE:=proc(Re1,Re2,Ak) local A,k,B,C,i,j,m,n,eq,J,sol,RE,RE1,RE2,coe,vars,rat,LCM,approach,dencoe; begin A:=op(Ak,0); k:=op(Ak,1); m:=OrderRE(Re1,A(k)); n:=OrderRE(Re2,A(k)); B:=genident();C:=genident(); coe:=[genident() $ i=1..n+m+1]; vars:=[B(k+i-1) $ i=1..m, C(k+i-1) $ i=1..n]; J:=max(m,n); RE1:=subs(lhs(Re1)-rhs(Re1),A=B); RE2:=subs(lhs(Re2)-rhs(Re2),A=C); approach:=0; for i from 1 to J-1 do rat:=ReduceRE(B(k+i-1),B(k),RE1)+ ReduceRE(C(k+i-1),C(k),RE2); dencoe[i]:=denom(rat); approach:=approach+coe[i]*numer(rat); end_for; for i from J to n+m+1 do rat:=ReduceRE(B(k+i-1),B(k),RE1)+ ReduceRE(C(k+i-1),C(k),RE2); dencoe[i]:=denom(rat); approach:=approach+coe[i]*numer(rat); eq:=numer(normal(approach)); eq:=collect(eq,vars,factor); eq:={coeff(eq,vars)}; sol:=solve(eq,{coe[j] $ j=1..i},IgnoreSpecialCases); if subs({coe[j] $ j=1..i},op(sol))={0} then else LCM:=lcm(dencoe[j] $ j=1..i); RE:=_plus(coe[j]*dencoe[j]/LCM*A(k+j-1) $ j=1..i); RE:=subs(RE,op(sol)); RE:=subs(RE,map(indets(RE) minus _union(indets(Re1),indets(Re2)),z->z=1)); RE:=numer(normal(RE)); RE:=collect(RE,[A(k+j-1) $ j=1..i],factor)=0; return(RE); end_if; end_for; end_proc: /*************************************************************************************************/ /* ProductRE determines a recurrence equation out of two recurrence equations Re1 and Re2, which is satisfied by the product of the solutions of Re1 and Re2 */ /*ProductRE:=proc(Re1,Re2,Ak) local A,k,B,C,DELTAf,i,j,l,m,n,eq,J,sol,RE,RE1,RE2,rules,productRE,vars,redRE; begin A:=op(Ak,0); k:=op(Ak,1); m:=OrderRE(Re1,A(k)); n:=OrderRE(Re2,A(k)); B:=genident();C:=genident(); productRE:=[genident() $ i=1..n*m+1]; vars:=[genident() $ i=1..n*m]; rules:=B(k+m-i)*C(k+n-j)=vars[m*(j-1)+i] $ i=1..m $ j=1..n; J:=max(m,n); RE1:=subs(lhs(Re1)-rhs(Re1),A=B); RE2:=subs(lhs(Re2)-rhs(Re2),A=C); eq:={}; redRE:=B(k)*C(k); for i from 1 to J-1 do eq:={op(eq),subsex(numer(normal(productRE[i]-expand(redRE))),rules)}; redRE:=DELTA(redRE,k); redRE:=ReduceRE(redRE,B(k),RE1); redRE:=ReduceRE(redRE,C(k),RE2); end_for; for i from J to n*m+1 do eq:={op(eq),subsex(numer(normal(productRE[i]-expand(redRE))),rules)}; sol:=groebner::eliminate(eq,vars); if nops(sol)<>0 then RE:=collect(op(sol),productRE,factor); DELTAf:=fp::fixargs(DELTA,1,k): RE:=subs(RE,productRE[l]=(DELTAf@@(l-1))(A(k)) $ l=1..i); RE:=collect(RE,[A(k+l-1) $ l=1..i],factor)=0; return(RE); end_if; redRE:=DELTA(redRE,k); redRE:=ReduceRE(redRE,B(k),RE1); redRE:=ReduceRE(redRE,C(k),RE2); end_for; end_proc:*/ ProductRE:=proc(Re1,Re2,Ak) local A,k,B,C,DELTAf,i,j,m,n,eq,J,sol,RE,RE1,RE2,rules,vars,redRE,LCM,approach,coe,dencoe,rat; begin A:=op(Ak,0); k:=op(Ak,1); m:=OrderRE(Re1,A(k)); n:=OrderRE(Re2,A(k)); B:=genident();C:=genident(); coe:=[genident() $ i=1..n*m+1]; vars:=[genident() $ i=1..n*m]; rules:=B(k+m-i)*C(k+n-j)=vars[m*(j-1)+i] $ i=1..m $ j=1..n; J:=max(m,n); RE1:=subs(lhs(Re1)-rhs(Re1),A=B); RE2:=subs(lhs(Re2)-rhs(Re2),A=C); approach:=0; redRE:=B(k)*C(k); for i from 1 to J-1 do rat:=normal(expand(redRE)); dencoe[i]:=denom(rat); approach:=approach+subsex(coe[i]*numer(rat),rules); redRE:=DELTA(redRE,k); redRE:=ReduceRE(redRE,B(k),RE1); redRE:=ReduceRE(redRE,C(k),RE2); end_for; for i from J to n*m+1 do rat:=normal(expand(redRE)); dencoe[i]:=denom(rat); approach:=approach+subsex(coe[i]*numer(rat),rules); eq:=numer(normal(approach)); eq:=collect(eq,vars,factor); eq:={coeff(eq,vars)}; sol:=solve(eq,{coe[j] $ j=1..i},IgnoreSpecialCases); if subs({coe[j] $ j=1..i},op(sol))={0} then else LCM:=lcm(dencoe[j] $ j=1..i); DELTAf:=fp::fixargs(DELTA,1,k): RE:=_plus(coe[j]*dencoe[j]/LCM*(DELTAf@@(j-1))(A(k)) $ j=1..i); RE:=subs(RE,op(sol)); RE:=subs(RE,map(indets(RE) minus _union(indets(Re1),indets(Re2)),z->z=1)); RE:=numer(normal(RE)); RE:=collect(RE,[A(k+j-1) $ j=1..i],factor)=0; return(RE); end_if; redRE:=DELTA(redRE,k); redRE:=ReduceRE(redRE,B(k),RE1); redRE:=ReduceRE(redRE,C(k),RE2); end_for; end_proc: /*************************************************************************************************/ /* DELTA represents the forward difference operator */ DELTA:=proc(a,k) local i,delta,deltalist,const; begin // Product if testtype(a,"_mult") then deltalist:={}; const:=1; for i from 1 to nops(a) do if has(op(a,i),k) then deltalist:={op(deltalist),op(a,i)} else const:=const*op(a,i) end_if; end_for; if nops(deltalist)<2 then return(const*DELTA(deltalist[1],k)) else delta:=deltalist[1]; for i from 2 to nops(deltalist) do delta:=DELTA(deltalist[i],k)*subs(delta,k=k+1)+ DELTA(delta,k)*deltalist[i] // product rule end_for; end_if; return(const*delta); end_if; // Sum if testtype(a,"_plus") then delta:=0; for i from 1 to nops(a) do delta:=delta+DELTA(op(a,i),k) // linearity end_for; return(delta); end_if; // Power if testtype(a,"_power") and testtype(op(a,2),Type::PosInt) then delta:=op(a,1); for i from 2 to nops(a) do delta:=DELTA(op(a,i),k)*subs(delta,k=k+1)+ DELTA(delta,k)*op(a,i) // product rule end_for; return(delta); end_if; subs(a,k=k+1)-a end_proc: /*************************************************************************************************/ /* OrderRE determines the order of a recurrence equation Re in AK */ OrderRE:=proc(Re,Ak) local A,k,RE,ind; begin A:=op(Ak,0); k:=op(Ak,1); if testtype(Re,"_equal") then RE:=lhs(Re)-rhs(Re); else RE:=Re end_if; ind:=map(select(indets(RE,RatExpr),has,A),x->(op(x,1)-k)); if nops(ind)>1 then max(ind) //-min(ind) /* to determine "real" order: unquote */ else op(ind) end_if; end_proc: /*************************************************************************************************/ /* ContentRE determines the content of a recurrence equation Re in Ak */ ContentRE:=proc(Re,Ak) local A,k,RE,i,content; begin A:=op(Ak,0); k:=op(Ak,1); if testtype(Re,"_equal") then RE:=lhs(Re)-rhs(Re); else RE:=Re end_if; RE:=factor(RE); content:=op(RE,1); for i from 2 to nops(RE) step 2 do if not(has(op(RE,i),A)) then content:=content*op(RE,i)^op(RE,i+1); end_if; end_for; content; end_proc: /*************************************************************************************************/ /* SolveRE solves a holonomic recurrence equation Re in Ak of hypergeometric type */ SolveRE:=proc(Re,Ak) local A,k,i,j,mn,mx,start,temp,RE,RE2,arg1,arg2,sol,rat,num,den,upp,low,minupp,minlow,coe0,coe1,x,order,shift,piecewisesol,content; begin A:=op(Ak,0); k:=op(Ak,1); RE:=lhs(Re)-rhs(Re); if testtype(RE,"_plus") and nops(RE)=2 then if nops(op(RE,1))=1 then arg1:=op(op(RE,1),1)-k; else arg1:=op(select(op(RE,1),has,A),1)-k; end_if; if nops(op(RE,2))=1 then arg2:=op(op(RE,2),1)-k; else arg2:=op(select(op(RE,2),has,A),1)-k; end_if; mx:=max(arg1,arg2); mn:=min(arg1,arg2); order:=mx-mn; RE:=subs(RE,k=k-mn); // call SolveRE recursively if order>1 then temp:=genident(); RE:=subs(RE,A(k+order)=A(temp+1),A(k)=A(temp)); piecewisesol:=[]; for i from 0 to order-1 do RE2:=RE; RE2:=subs(RE2,k=order*k+i); RE2:=subs(RE2,temp=k); sol:=SolveRE(RE2=0,A(k)); if op(sol,0)=A then start:=op(sol,1); else start:=op(select(sol,has,A),1); end_if; sol:=subs(sol,A(start)=A(order*start+i)); sol:=subs(sol,k=(k-i)/order); piecewisesol:=[op(piecewisesol),[_mod(k,order)=i,sol]] end_for; return(piecewise(op(piecewisesol))) end_if; sol:=normal(op(solve(RE,A(k+order),IgnoreSpecialCases))/A(k)); rat:=simplify(expand(sol)); num:=factor(numer(rat)); den:=factor(denom(rat)); x:=1; // determine upper pochhammer symbols upp:=[];minupp:=infinity; if nops(num)>2 then x:=x*op(num,1) else x:=x*num end_if; for i from 2 to nops(num) step 2 do coe1:=coeff(op(num,i),k,1); coe0:=coeff(op(num,i),k,0); if coe1<>0 then x:=x*coe1^op(num,i+1); if testtype(coe0/coe1,Type::Integer) then minupp:=min(minupp,coe0/coe1) end_if; for j from 1 to op(num,i+1) do upp:=[op(upp),coe0/coe1] end_for; else x:=x*op(num,i)^op(num,i+1) end_if; end_for; // determine lower pochhammer symbols low:=[];minlow:=infinity; if nops(den)>2 then x:=x/op(den,1) else x:=x/den end_if; for i from 2 to nops(den) step 2 do coe1:=coeff(op(den,i),k,1); coe0:=coeff(op(den,i),k,0); if coe1<>0 then x:=x/coe1^op(den,i+1); if testtype(coe0/coe1,Type::Integer) then minlow:=min(minlow,coe0/coe1) end_if; for j from 1 to op(den,i+1) do low:=[op(low),coe0/coe1] end_for; else x:=x/op(den,i)^op(den,i+1) end_if; end_for; // determine shift content:=ContentRE(RE,A(k)); shift:=min(minlow,minupp,op(map(solve(content,k),z->-z))); if testtype(shift,Type::NegInt) or shift=0 then shift:=shift-1; upp:=map(upp,y->y-shift); low:=map(low,y->y-shift); else shift:=0 end_if; return(subs(A(-shift)*x^k* _mult(op(map(upp,y->pochhammer(y,k))))/ _mult(op(map(low,y->pochhammer(y,k)))),k=k+shift)) else error("Input is no recurrence equation of hypergeometric type") end_if; end_proc: /*************************************************************************************************/ /* FPS determines a representation as formal power series (with summation index k) for a function f in x */ FPS:=proc(f,xeq,k,simplified=TRUE) local DE,RE,F,A,x,x0,n,m,i,coepol,ratzero,lcm,fps,sol,shift,st; begin F:=genident("F"); A:=genident("A"); if testtype(xeq,"_equal") then x:=lhs(xeq); x0:=rhs(xeq); else x:=xeq; x0:=0; end_if; st:=time(); DE:=HolonomicDE(f,F(x)); userinfo(3,"holonomic differential equation: \n".expr2text(DE)); userinfo(2,"HolonomicDE: ".expr2text(time()-st)." ms"); /* st:=time(); DE:=SubsDE(DE,F(x),x+x0); userinfo(3,"holonomic differential equation: \n".expr2text(DE)); userinfo(2,"SubsDE: ".expr2text(time()-st)." ms"); */ st:=time(); RE:=DEtoRE(DE,F(x),A(k)); userinfo(3,"holonomic recurrence equation: \n".expr2text(RE)); userinfo(2,"DEtoRE: ".expr2text(time()-st)." ms"); n:=OrderRE(RE,A(k)); if n=1 then coepol:=coeff(lhs(RE),A(k+1)); ratzero:=solve(coepol,k); ratzero:=ratzero minus select(ratzero,testtype,Type::Integer); if nops(ratzero)>0 then // call FPS recursively lcm:=ilcm(op(map(ratzero,y->denom(y)))); fps:=FPS(subsex(f,x=x^lcm,(x^lcm)^(1/lcm)=x),x,k); return(subs(fps,x=x^(1/lcm),(x^(1/lcm))^lcm=x)); else st:=time(); sol:=SolveRE(RE,A(k)); userinfo(3,"solution of holonomic recurrence equation: \n".expr2text(sol)); userinfo(2,"SolveRE: ".expr2text(time()-st)." ms"); if nops(sol)=1 then shift:=op(sol,1); else shift:=op(select(sol,has,A),1); end_if; sol:=subs(sol*x^k,k=k+shift); sol:=subs(sol,A(shift)=limit(diff(f,x$shift),x=0)/shift!); st:=time(); if simplified then sol:=Simplify(sol) end_if; userinfo(2,"Simplify: ".expr2text(time()-st)." ms"); return(SUM(sol,k=0..infinity)) end_if; end_if; st:=time(); sol:=SolveRE(RE,A(k)); userinfo(3,"solution of holonomic recurrence equation: \n".expr2text(sol)); userinfo(2,"SolveRE: ".expr2text(time()-st)." ms"); m:=nops(sol); shift:=[]; for i from 1 to m do if nops(piecewise::expression(sol,i))=1 then shift:=[op(shift),op(piecewise::expression(sol,i),1)] else shift:=[op(shift),op(select(piecewise::expression(sol,i),has,A),1)] end_if; end_for; sol:=[subs(piecewise::expression(sol,i)*x^k,k=m*k+shift[i]) $ i=1..m]; sol:=subs(sol,{A(shift[i])=limit(diff(f,x$shift[i]),x=0)/shift[i]! $ i=1..m}); st:=time(); if simplified then sol:=map(sol,y->Simplify(y)) end_if; userinfo(2,"Simplify: ".expr2text(time()-st)." ms"); subs(_plus(SUM(sol[i],k=0..infinity) $ i=1..m),SUM(0,k=0..infinity)=0) end_proc: /*************************************************************************************************/ /* TPS determines a representation as truncated power series for a function f in x */ TPS:=proc(f,x,ord) local DE,RE,F,A,k,i,rules; begin F:=genident(); A:=genident(); k:=genident(); DE:=HolonomicDE(f,F(x)); RE:=DEtoRE(DE,F(x),A(k)); n:=OrderRE(RE,A(k)); RE:=op(solve(RE,A(k+n),IgnoreSpecialCases)); rules:={A(i)=limit(diff(f,x$i),x=0)/i! $ i=0..n-1}; // initial values for i from n to ord do rules:={op(rules),A(i)=subs(RE,k=i-n,rules)} end_for; subs(_plus(A(i)*x^i $ i=0..ord),rules) end_proc: /*************************************************************************************************/ print(NoNL,"Package 'Formal Power Series 1.0' , MuPAD 3.1.1\n"); print(NoNL,"Torsten Sprenger, University of Kassel , 2006\n");