Artifact fe07f96e3c5dae377f0102745407cc0d67a04bd9f09c75c3cc592b55a31322ac:
- Executable file
r36/src/fps.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 36329) [annotate] [blame] [check-ins using] [more...]
module fps; % a package for the algorithmic calculation % of Formal Power Series % Author : Wolfram Koepf, ZIB Berlin % REDUCE version by: Winfried Neun, ZIB Berlin. % (September 1994) % |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| % % % % Please report bugs to Winfried Neun, % % Konrad-Zuse-Zentrum % % fuer Informationstechnik Berlin, % % Heilbronner Str. 10 % % 10711 Berlin - Wilmersdorf % % Federal Republic of Germany % % or by email, neun@sc.ZIB-Berlin.de % % % % |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| % create!-package ('(fps simplede substexp linrec hgrsolve constre ratalgo), '(contrib specfn)); load_package misc; % For limits. load_package factor,specfn; fluid '(ps!:order!-limit); lisp(ps!:order!-limit := 30); algebraic << factor factorial >>; endmodule; module simplede; fluid '(fps_search_depth !*protfg ps!:order!-limit); global '(inconsistent!*); share fps_search_depth; fps_search_depth := 5; %the default switch tracefps; algebraic << operator ba; operator infsum ; array dff(50) >>; put('simplede,'psopfn,'simpledeeval); symbolic procedure simpledeeval(u); begin scalar res,usevar; if length u = 2 then << usevar := 'y; res := int_simplede(car u,cadr u); if eq(res,-1) then return simpledeexit(car u,cadr u,'y); >> else if length u = 3 then << usevar := caddr u; res := int_simplede(car u,cadr u); if eq(res,-1) then return simpledeexit(car u,cadr u,usevar); >> else rederr("Wrong number of Arguments for simplede"); res := sublis('((oddexpt . expt)(ba . a)(nn . k)),res); return if reval usevar = usevar then sublis(list('ff . usevar),res) else sublis(list('ff . intern gensym()),res); end; algebraic procedure int_simplede(f,x); begin scalar cap_a,degree0fde,cap_f,j,cap_j,nnn,s,ind,deq,eqq,reqq, ak,terms,list1,list2,nmax,cap_m,cap_r,ii,m,leadcoeff,m0, len,cap_s,result,parameters,solved,!*allfac,!*protfg; !*protfg := t; nmax :=fps_search_depth; clear a; operator a; off allfac; depend ff,x; dff(0) := f; % start search for a simple DE for degreeofde:=1:nmax do << dff(degreeofde) := df(dff(degreeofde-1),x); eqq := dff(degreeofde) + for j:=0:(degreeofde-1) sum a(j) * dff(j); eqq := recursionsimplify(eqq); eqq := num eqq; terms := {}; list1 := converttolist (eqq,degreeofde+1); while list1 neq {} do << list2 := {}; j := fastpart(list1 ,1); cap_j := j; len := fastlength list1; for i:=2:len do if type_ratpoly(j/fastpart(list1,i),x) then cap_j := cap_j + fastpart(list1,i) else list2 := fastpart(list1,i) . list2; terms := cap_j . terms; list1 := reverse list2; >>; ind := for j:=0:degreeofde-1 collect a(j); s := savesolve(terms,ind); if s = {} then nil else << if symbolic !*tracefps then write "Solution: ",s; result := degreeofde; nmax := 0 >>; >>; degreeofde := result; if nmax = 0 then << if symbolic !*tracefps then write " successful search for DE">> else return -1; for each ss in first s do << ss := sub(a(degreeofde)=1,ss); setze(lhs ss,rhs ss)>>; % setting up the Differential equation on factor; deq := df(ff,x,degreeofde) + for j:=0:(degreeofde-1) sum a(j)*df(ff,x,j); off factor; deq := num deq; return deq; end; put('fps,'psopfn,'fpseval); symbolic procedure fpseval(u); begin scalar gens,res,!*factor,!*precise; if length u = 2 then << res := psalg(car u,cadr u); if eq(res,-1) then return fpsexit(car u,cadr u,0); return sublis('((oddexpt . expt)(ba . a) (nn . k)),res) >> else if length u = 3 then << gens := gensym(); res := psalg(sublis(list( cadr u . gens),car u),gens); if eq(res,-1) then return fpsexit(car u,cadr u,caddr u); res := sublis('((oddexpt . expt)(ba . a) (nn . k)),res); res := subf(caadr res, list list(gens,'plus,cadr u, list('minus,caddr u))); return mk!*sq res; >> else rederr("Wrong number of Arguments for FPS"); end; algebraic procedure asymptpowerseries (f,x); sub(x=1/x,fps(sub(x=1/x,f),x)); symbolic procedure fpsexit(a,b,z); << erfg!* := nil; list ('fps,a,b, z) >>$ symbolic procedure simpledeexit(a,b,z); << erfg!* := nil; list ('simplede,a,b, z) >>$ algebraic procedure psalg(f,x); begin scalar cap_a,degree0fde,cap_f,j,cap_j,nnn,s,ind,deq,eqq,reqq, ak,terms,list1,list2,nmax,cap_m,cap_r,ii,m,leadcoeff,m0, len,cap_s,result,parameters,solved,!*allfac,!*protfg; f := recursionsimplify f; !*protfg := t; nmax :=fps_search_depth; clear a; operator a; off allfac; depend ff,x; dff(0) := f; % special cases if polynomq(f,x) then return f; if type_ratpoly(f,x) then return ratalgo(f,x); % start search for a simple DE clearrules special!*pochhammer!*rules; clearrules spec_factorial; clearrules spec_pochhammer; for degreeofde:=1:nmax do << dff(degreeofde) := df(dff(degreeofde-1),x); eqq := dff(degreeofde) + for j:=0:(degreeofde-1) sum a(j) * dff(j); eqq := recursionsimplify(eqq); eqq := num eqq; terms := {}; list1 := converttolist (eqq,degreeofde+1); while list1 neq {} do << list2 := {}; j := fastpart(list1,1); cap_j := j; len := fastlength list1; for i:=2:len do if type_ratpoly(j/fastpart(list1,i),x) then cap_j:= cap_j + fastpart(list1,i) else list2 := fastpart(list1,i) . list2; terms := cap_j . terms; list1 := reverse list2; >>; ind := for j:=0:degreeofde-1 collect a(j); s := savesolve(terms,ind); if s = {} then nil else << if symbolic !*tracefps then write "Solution: ",s; result := degreeofde; nmax := 0 >>; >>; degreeofde := result; if nmax = 0 then << if symbolic !*tracefps then write " successful search for DE">> else return -1; for each ss in first s do << ss := sub(a(degreeofde)=1,ss); setze(lhs ss,rhs ss)>>; % setting up the Differential equation on factor; deq := df(ff,x,degreeofde) + for j:=0:(degreeofde-1) sum a(j)*df(ff,x,j); off factor; deq := num deq; if symbolic !*tracefps then write("Differential equation is: ", deq); % transforming into Recurrence equation factor ba; let subst_rules; req := pssubst(deq,x,ba,nn); clearrules subst_rules; if symbolic !*tracefps then write("Recurrence equation is :",req); ind := {}; for ii:=-50 : 50 do if not(coeffn(req,ba(nn+ii),1) =0) then ind := ii . ind; cap_m := first ind; if symbolic !*tracefps then write(" M, ind, parameters : ",cap_m,",",ind,",", parameters); leadcoeff := num coeffn(req,ba(nn+cap_m),1); nnn := fastlength ind; let special!*pochhammer!*rules; let spec_factorial; let spec_pochhammer; result := 0; if (nnn = 1) then << % functions with finite representation if symbolic !*tracefps then write "fps with finite number of non-zero coefficients"; cap_r := sub(nn=nn+(1-cap_m),- reduct(req,ba(nn+cap_m))) /(sub(nn=nn+(1-cap_m),lcof(req,ba(nn+cap_m))) * ba(nn)); leadcoeff:= sub(nn=nn+(1-cap_m),leadcoeff); result := constantre(cap_r,leadcoeff,0,nn,x); if result = failed then result :=0; >>; % try hypergeometric case if (nnn = 2) then << m := abs(first ind - second ind); cap_r := sub(nn=nn+(m-cap_m),- reduct(req,ba(nn+cap_m))) /(sub(nn=nn+(m-cap_m),lcof(req,ba(nn+cap_m))) * ba(nn)); leadcoeff:= sub(nn=nn+(m-cap_m),leadcoeff); result := hypergeomre(m,cap_r,leadcoeff,0,nn,x) >>; if result =0 then << % test for constant coefficients terms := for j:=0:(degreeofde-1) join if freeof(a(j),x) then {} else {t}; if terms = {} then << req := ba(k+degreeofde) + for j:=0:(degreeofde-1) sum ba(k+j)*a(j); if symbolic !*tracefps then << write("DE has constant coefficients"); write("DE = ",deq); write("RE = ",req); >>; s := 0; iii := 0; while freeof(req,ba(k + iii)) do << s := s + limit(dff(iii),x,0) * x^iii; iii := iii + 1 >>; m0 := iii; if symbolic !*tracefps then write "i was found : ",iii; if m0 <= degreeofde-1 then << s := solve_lin_rec(req,for i:=m0:(degreeofde-1) collect ba(i) = limit(dff(i),x,0)); if symbolic !*tracefps then write("solution : ",s); s:=sub(n=nn,s); result := infsum(s/(factorial nn) * x^nn,nn,0,infinity) >> else result := s; >>; >>; if result = 0 or not(freeof(result,failed)) then return (-1); lisp (erfg!* := nil); result:= result; let hgspec_pochhammer; result:= result; clearrules hgspec_pochhammer; result := verbessere (result,nil); return result; end; flag ('(verbessere), 'opfn); symbolic procedure verbessere (x,uu); << if eqcar (x,'plus) then 'plus . foreach xx in cdr x collect verbessere(xx , nil) else if not (eqcar (x,'infsum)) then x else << if eqcar (x,'infsum) and eqcar(cadr x,'quotient) then x := list('infsum ,list('quotient, simplify_expt cadr cadr x,simplify_expt caddr cadr x)); uu := cadr x; if eqcar (x,'infsum) and eqcar(cadr x,'quotient) then << uu := int_simplify_factorial auxcopy cadr x >>; list('infsum , uu,'nn,0,'infinity)>> >>; symbolic procedure zerlege u; if fixp u and u>0 and (u<10000 or !*ifactor) then << u := zfactor u; for each j in u join for jj :=1:cdr j collect (car j) >> else list(u); symbolic procedure simplify_expt u; begin scalar uu,exptlist,nonexptlist,asso,numb,expo; uu := u; if eqcar(u,'times) then u := cdr u; while u do << if pairp car u and (eq (caar u,'expt) or eq (caar u,'sqrt)) then << if numberp cadar u then numb := zerlege (cadar u) else numb := list cadar u; expo := if eq (caar u,'sqrt) then '((quotient 1 2)) else cddar u; while numb do << if asso:= atsoc (car numb,exptlist) then exptlist := (car numb . list list('plus,car expo,cadr asso)) . delasc (car numb,exptlist) else exptlist := ((car numb) . expo) . exptlist ; numb := cdr numb; >>; >> else if and(idp car u,asso := atsoc (car u,exptlist)) then << exptlist := (car u . list list('plus,1,cadr asso)) . delasc (car u,exptlist) >> else nonexptlist := (car u) . nonexptlist; u := cdr u; >>; if null exptlist then return uu; for each x in exptlist do nonexptlist := ('oddexpt . x ) . nonexptlist; return (car uu) . nonexptlist; end; fluid ('(rsolve!*!*)); algebraic procedure hypergeomre(m,cap_r,leadcoeff,dffpointer,k,x); % solve the hypergeometric Recurrence Equation % % a(k+m) = cap_R(k) * a(k) % % where leadcoeff is the leading coefficient of the RE % and DF is a table where DF(dffpointer+i) = df(f,x,i) begin scalar denr,fract,ii,m0,m1,c0,ck,s,c,df2,q,r2,lterm,nn, s0, leadcoeff2; denr := solve(leadcoeff,k); m0 := {}; foreach xx in denr do if type_rational rhs xx then m0 := ((rhs xx)+1) . m0; if not(m0 = {}) then m0 := max(m0) else m0 := 0; if symbolic !*tracefps then << write "RE is of hypergeometric type"; write "Symmetry number mm := ",m; write "RE: for all k >= ",m0,": a (k + ",m,") = " ,cap_r * a(k); write "leadcoeff := ",leadcoeff; >>; fract := {}; foreach xx in denr do if type_fraction(rhs xx) then fract := den(rhs xx) . fract; if not(fract = {}) then << q := first fract; dff(dffpointer + 10) := sub(x=x^q,dff(dffpointer)); if symbolic !*tracefps then << write "RE modified to nn= ",k/q; write "=> f := ",dff(dffpointer + 10)>>; s := hypergeomre(q*m,sub(k=k/q,cap_r), sub(k=k/q,leadcoeff),dffpointer + 10,k,x); return sub(x=x^(1/q),s); >>; if m0 < 0 then << nn:= -m0 + remainder(-m0,m); dff(dffpointer + 10) := df2 := x^nn * dff(dffpointer); if symbolic !*tracefps then << write "working with ",x^nn,"*f"; write "=> f :=" ,df2 >>; s := hypergeomre(m,sub(k=k-nn,cap_r), sub(k=k-nn,leadcoeff), dffpointer + 10,k,x); return update_coeff(s,x,-nn) >>; if m0 > 0 then << m1 := {}; foreach xx in denr do if type_rational rhs xx then m1 := append(list (rhs xx +1),m1); m1 := min m1; if m1 > 0 then << dff(dffpointer + 10) := df2 := x^(-m1)*dff(dffpointer); if symbolic !*tracefps then << write"a(k) = 0 for k < ",m1; write "working with ",x^(-m1),"*f"; write "=> f :=" ,df2 >>; s := hypergeomre(m,sub(k=k+m1,cap_r), sub(k=k+m1,leadcoeff), dffpointer + 10,k,x); return update_coeff(s,x,m1); >> >>; % logarithmic singularity Baustelle if lisp pairp errorset!*(list ('simptaylor,mkquote list( mkquote list('dff,dffpointer), mkquote x,0,1)),nil) then << lterm := num taylortostandard(taylor(dff(dffpointer),x,0,1)); nn := 0; if lisp(if member('(log x) ,kernels !*q2f simp lterm) then t else nil) % Comments? then << dff(dffpointer + 10):=dff(dffpointer) -lterm; if symbolic !*tracefps then write "=> f :=",dff(dffpointer + 10); s := hypergeomre(m, r, leadcoeff*(k-nn), dffpointer + 10,k,x); return(lterm+s); >>; >>; s := 0; s0 := 0; for i:=0:(m0+m-1) do << if i > 0 then dff(dffpointer + i) := df(dff(dffpointer + i-1),x); c0 := limit(dff(dffpointer + i),x,0); if (lisp listp reval c0 and fastpart(c0,0) = limit) then << if symbolic !*tracefps then write "Could not find the limit of: " ,dff(dffpointer + i),",",x,",",0; rederr("Problem using limit operator") >> else << c0 := c0/factorial (i); if symbolic !*tracefps then write " a(",i,") = ",c0; if not (c0 =0) then << s0 := s0+c0*x^i; if i < m0 then s := s + c0 * x^i % single terms else << ck := hypergeomrsolve(sub(k=m*k+i,cap_r),k,c0); if symbolic !*tracefps then write " ck = ",ck; c :=1; ck := ck/c; let hgspec_pochhammer; ck := ck; clearrules hgspec_pochhammer; if ck = 0 then s := s + c0*x^i else if rsolve!*!* = finite then s := s + c*sum(ck*x^(m*k+i), k) else s := s + c * infsum(ck*x^(m*k+i)) ; if symbolic !*tracefps then write " S = ",s; >> >> >> >>; return (s); end; algebraic << let infsum(0) = 0>>; % some compatibility functions for Maple sources. symbolic flag('(savesolve type_fraction type_rational),'opfn); algebraic procedure converttolist (express,len); << len := fastlength express; for i:=1:len collect fastpart(express , i)>>; symbolic procedure type_fraction (num); (if pairp num1 and fixp car num1 and fixp cdr num1 and not onep cdr num1 then num else nil) where num1 := simp num; symbolic procedure type_rational(num); (if pairp num1 and (fixp car num1 or null car num1) and fixp cdr num1 then t else nil) where num1 := simp num; algebraic procedure type_ratpoly(exprn,var); if (polynomq (den exprn,var) and polynomq (num exprn,var)) then t else nil; symbolic procedure savesolve (x,y); begin scalar !*cramer; on cramer; % this is a temporary fix for solve !! % check with fps(sin(x)^2 * cos(x)^2,x); !! return << switch solveinconsistent; on solveinconsistent; inconsistent!* := nil; if pairp (x := errorset!*(list ('solveeval,mkquote list(x,y)),nil)) and not inconsistent!* then << x :=car x; if x = '(list) then x else if eqcar(cadr x,'equal) then % one element solution list('list,x) else x>> else list('list) >>; end; algebraic procedure setze(x,y); let x = y; symbolic procedure polynomq (x,var); if not fixp denr simp x then nil else begin scalar kerns,kern; kerns := kernels !*q2f simp x; aa: if null kerns then return t; kern := first kerns; kerns := cdr kerns; if not(eq (kern, var)) and depends(kern,var) then return nil else go aa; end; flag('(polynomq),'opfn); flag ('(polynomq type_ratpoly),'boolean); algebraic << operator update_coeff; update_coeff_rules := { update_coeff (~a + ~b,~x,~m) => update_coeff(a,~x,~m) + update_coeff(b,~x,~m), update_coeff (~c * ~a,~x,~m) => c * update_coeff(a,~x,~m) when freeof(c,x), update_coeff ( - ~a,~x,~m) => - update_coeff(a,~x,~m), update_coeff (~a/~c,~x,~m) => update_coeff(a,~x,~m) /c when freeof(c,x) and c neq 1, update_coeff (~x,~x,~m) => x^(m + 1), update_coeff (~c,~x,~m) => c * x^m when freeof(c,x), update_coeff (infsum(~xx),~x,~m) => infsum update_coeff(xx,x,m), update_coeff (~x^~j*~xx,~x,~m) => x^(j + m + 1)when x = xx, update_coeff (~x^~j*~xx^~jj,~x,~m) => x^(j + jj + m) when x = xx, update_coeff (~x^~j,~x,~m) => x^(j + m)}$ let update_coeff_rules >>$ endmodule; module substexp; algebraic; depend ff,x; depend f,x; depend f,!~k; operator pssubst; operator a,pproduct; subst_rules := { pssubst(- ~g,~x,~a,~n) => -pssubst(g,x,a,n), pssubst(~g+~h,~x,~a,~n) => pssubst(g,x,a,n) + pssubst(h,x,a,n), pssubst(~c*~g,~x,~a,~n) => c*pssubst(g,x,a,n) when freeof(c,x) and freeof(c,g), pssubst(df(~f,~x,~k),~x,~a,~n) => pochhammer(n+1,k) * a(n+k), pssubst(df(~f,~x),~x,~a,~n) => (n + 1)* a(n + 1), pssubst(~x^~j * df(~f,~x),~x,~a,~n) => pochhammer(n+1-j,1)*a(n+1-j), pssubst(~x^~j * df(~f,~x,~k),~x,~a,~n)=> pochhammer(n+1-j,k)*a(n+k-j), pssubst(ff*~x^~j,~x,~a,~n) => a(n-j), pssubst(ff*~x,~x,~a,~n) => a(n-1), pssubst(df(~f,~x) *x,~x,~a,~n) => pochhammer(n,1)*a(n), pssubst(df(~f,~x,~k) *x,~x,~a,~n) => pochhammer(n,k)*a(n+k-1), pssubst(f,~x,~a,~n) => a(n), pssubst(ff,~x,~a,~n) => a(n), pssubst(~c,~x,~a,~n) => 0 when freeof(c,x), pssubst(~x^~j,~x,~a,~n) => 0 when fixp j, pssubst(~x,~x,~a,~n) => 0 }; spec_pochhammer := { pochhammer(~a,~k)//pochhammer(~b,~k) => (a + k -1)/(a - 1) when (a - b)=1, pochhammer(~a,~k)//pochhammer(~b,~k) => (b - 1)/(b + k -1) when (b - a)=1, pochhammer(~z,~k) * pochhammer(~cz,~k) => prod((repart(z) + (j - 1))^2 + (impart(z))^2,j,1,k) when not(impart(z) = 0) and z = conj cz, pochhammer(~k,~n) => 1 when n=0, pochhammer(~k,~n) => pproduct (k,n) when fixp n, pproduct (~k,~ii) => 1 when ii =0, pproduct (~k,~ii) => (k + ii - 1) * pproduct (k,ii -1)}$ spec_factorial := { factorial (~n) // factorial (~n+1) => 1/(n+1), factorial (~n) * factorial (~n) // factorial (~n+1) => factorial (n)/(n+1), factorial (~n+1) // factorial (~n) => (n+1), factorial (~n+1) * factorial (~n+1) // factorial (~n) => (n+1) * factorial (~n+1), (~otto ^(~k)) * factorial (~n) // factorial (~n +1) => otto^k /(n+1), (~otto ^(~k)) * factorial (~n+1) // factorial (~n) => otto^k * (n+1), (~otto ^~k) * ~hugo * factorial (~n) // factorial (~n +1) => otto^k * hugo/(n+1), (~otto ^~k) * ~hugo * factorial (~n+1) // factorial (~n) => otto^k * hugo *(n+1)}$ endmodule; module linrec; % solves (simple in) homogenous linear recursion relation REC which % has to have the form var(n)... with initial conditions % IC in the form of list {var(1)=?,var(2)=...} % the method used - substitution s(n)=x^n % like e^(lambda*x) in linear homogenous ODE. % The inhomogenous term has to be a constant. % The following code has been developed by Richard Liska from Prague % for the "Barry Simon Tests for PC Magazine". % Some checks and generalizations included by W. Neun , ZIB Berlin % the switch trlinrec turns the verbose mode on. switch trlinrec; algebraic procedure rsolve (rec,ic); begin scalar mvar1,inde,constpart,lowestind,modic,modrec,indis; mvar1 := part(mainvar rec ,0); inde := mainvar part(mainvar rec,1); indis := the_indices(foreach kk in ic collect lhs kk,mvar1); highestind := if indis neq {} then max(indis) else nil; nindis := the_indices(rec,mvar1); % applying rule : rec where r(~n) => 0 constpart := lisp aeval list('whereexp, list('list,list('replaceby,list(mvar1,list('!~,'n)),0)), rec); if not freeof(constpart,inde) then rederr ("Cant solve recurrence equations with non-constant coefficients"); if eq(constpart , 0) then return solve_lin_rec(rec,ic) else << modrec := sub (inde = inde +1,rec); modrec := modrec -rec; if (highestind) then << % Propagate the recursion to get additional start value modic := sub(inde=inde -max(nindis)+ highestind+1,rec); for each aa in ic do modic := sub(aa,modic); modic := (first solve(modic,mainvar modic)) . ic >> else modic := ic; return solve_lin_rec(modrec,modic); >>; end; fluid '(linrecx!*,linrecvar!*); algebraic procedure solve_lin_rec(rec,ic); % solves homogenous linear recursion relation REC which % has to have the form var(n)... with initial conditions % IC in the form of list {var(1)=?,var(2)=...} % the method used - substitution s(n)=x^n % like e^(lambda*x) in linear homogenous ODE. % Of course some checking should be added. (Done WN) begin scalar lrec,sol,msol,gsol,j,flagg,c,linrecvar!*,errflag,nsave; clear n; linrecvar!* := part (mainvar rec,0); linrecx!* := lisp gensym(); c:= lisp mkquote gensym(); %this is the dirty part. WN operator c; if(part(rec,0) eq linrecvar!*) and arglength(ic)=1 and part(mainvar lhs first ic,0) = linrecvar!* then return rhs (first ic); if(part(rec,0) neq plus) then return rederr "Cant solve recurrence equations with non-constant coefficients"; lrec := arglength rec; lrec := part(rec,lrec); for all n let linrecvar!*(n) = linrecx!*^n; lrec := lrec; rec:=rec /lrec; for all n clear linrecvar!* (n); rec:=num rec; for each j in coeff(rec,linrecx!*) do if (not freeof(j,part(part (part rec,1),1))) then errflag := 17; %??? if (errflag = 17) then return rederr "Cant solve recurrence equations with non-constant coefficients"; j:=1; for each a in solve(rec,linrecx!*) do <<a:=rhs a^n; gsol:=gsol+c(j)*a; j:=j+1; msol:=first multiplicities!*; multiplicities!*:=rest multiplicities!*; while msol>1 do <<a:=n*a; gsol:=gsol+c(j)*a; j:=j+1; msol:=msol-1 >> >>; if lisp !*trlinrec then write "General solution: ",linrecvar!* ,"(N) := ",gsol; if ic = {} then sol := {} else sol:=solve(for each a in ic collect sub(n=part(lhs a,1),gsol)=rhs a, for i:=1:arglength ic collect c(i)); % If some c(i) remains it can be arbitrary complex; sol := lisp subla('((equal . replaceby)),sol); sol := lisp subla('((equal . replaceby)),sol); let sol; gsol:=gsol; clearrules sol; let moivre_expt; gsol:=gsol; clearrules moivre_expt; for i:=1:j do if coeff(gsol,lisp list(c, i)) = list(gsol) then nil else gsol:= sub(lisp list(c, i) = lisp caaar makearbcomplex(),gsol); return gsol end; % (1 + i)**n => %applying Moivre's formula for complex numbers % N/2 N*PI N*PI % 2 *(COS(------) + SIN(------)*I) % 4 4 algebraic (moivre_expt := { (~z)^(~k) => moivre(z,k) when not freeof(z,i)}); algebraic procedure moivre(z,k); begin scalar rho,phi; % what ( will happen rho := sqrt( (repart z)^2 + (impart z)^2); if repart z = 0 then phi := pi/2 else phi := atan((impart z)/(repart z)); return rho^k *(cos(k*phi) + i * sin (k*phi)); end; algebraic procedure the_indices(ex,mvar); if part(ex,0) = list then for each kk in ex join the_indices(kk,mvar) else begin scalar eqq,l1,l2,kern; eqq := ex; lisp (kern := union (kernels !*q2f (numr simp eqq ./ 1), kernels !*q2f (denr simp eqq ./ 1))); l1 := 'list . lisp foreach k in kern join if atom k then nil else if eqcar(k,mvar) then list cadr k else nil; return l1; end; endmodule; module hypergeomrsolve; fluid '(!*tracefps); algebraic procedure hypergeomrsolve (r,k,a0); % solves the recurrence equation % % a(k+1) = r(k) * a(k), a(0) = a0 begin scalar re,nnn,ddd,c,p,q,ak,sols,ii; p := {}; q := {}; c := 1; re := r * (k + 1); nnn := factorize num re; ddd := factorize den re; foreach nn in nnn do if freeof (nn,k) then c := c * nn else if deg(nn,k) =1 then << c:= c*coeffn(nn,k,1); p:= append (p,list(coeffn(nn,k,0)/coeffn(nn,k,1)))>> else if deg(nn,k) =2 then << c := c * lcof(nn,k); sols := solve(nn,k); for each s in sols do << for i:=1:first multiplicities!* do p:= (- rhs s) . p; multiplicities!* := rest multiplicities!*; >> >> else rederr(" hypergeomRsolve failed"); foreach dd in ddd do if freeof (dd,k) then c := c / dd else if deg(dd,k) =1 then << c:= c/coeffn(dd,k,1); q:= append (q,list(coeffn(dd,k,0)/coeffn(dd,k,1)))>> else if deg(dd,k) =2 then << c := c / lcof(dd,k); sols := solve(dd,k); for each s in sols do << for i:=1:first multiplicities!* do q:= (- rhs s) . q; multiplicities!* := rest multiplicities!*; >>; >> else rederr(" hypergeomRsolve failed"); rsolve := infinite; for each s in p do if fixp s and s < 0 then rsolve := finite; if symbolic !*tracefps then write "RSOLVE = ",rsolve; p := for each s in p product pochhammer(s,k); q := for each s in q product pochhammer(s,k); ak := a0 * (c^k) * p/(q * factorial k); % Do additional simplification here?? return ak; end; % A special ruleset for powerseries; nn has a special meaning here and % should be treated as integer algebraic << hgspec_pochhammer := { pochhammer(~kk,~nn) => 1 when nn=0, pochhammer(~kk,nn) => 0 when kk = 0, pochhammer(~kk,nn) => (-1)^nn * factorial(-kk)/factorial(-kk-nn) when fixp(kk) and kk <=0, pochhammer(~kk,nn) => factorial(kk+nn-1)/factorial(kk-1) when fixp kk, pochhammer(~kk,~w*nn) => factorial(kk+w*nn-1)/factorial(kk-1) when fixp kk} >>; endmodule; module constre; algebraic procedure constantre(cap_r,leadcoeff,dffpointer,k,x); % solve constant RE % % a(k+1) = cap_R(k) * a(k) % % where leadcoeff is the leading coefficient of the RE % and DF is a table where DF(dffpointer+i) = df(f,x,i) begin scalar denr,fract,ii,m0,m1,c0,ck,s,c,df2,q,r2,lterm,nn, s0, leadcoeff2; denr := solve(leadcoeff,k); m0 := {}; foreach xx in denr do if type_rational rhs xx then m0 := ((rhs xx)+1) . m0; if not(m0 = {}) then m0 := max(m0) else m0 := 0; if symbolic !*tracefps then << write ">>> m0 = ",m0; write "RE is constant"; write "RE: for all k >= ",m0,": a (k + 1) = " ,cap_r * a(k); write "leadcoeff := ",leadcoeff; >>; fract := {}; foreach xx in denr do if type_fraction(rhs xx) then fract := den(rhs xx) . fract; if not(fract = {}) then << q := first fract; dff(dffpointer + 10) := sub(x=x^q,dff(dffpointer)); if symbolic !*tracefps then << write "RE modified to nn= ",k/q; write "=> f := ",dff(dffpointer + 10)>>; s := constantre(sub(k=k/q,cap_r), sub(k=k/q,leadcoeff),dffpointer + 10,k,x); return sub(x=x^(1/q),s); >>; if m0 < 0 then << nn:= -m0 + 1; dff(dffpointer + 10) := df2 := x^nn * dff(dffpointer); if symbolic !*tracefps then << write "working with ",x^nn,"*f"; write "=> f :=" ,df2 >>; s := constantre(sub(k=k-nn,cap_r),sub(k=k-nn,leadcoeff), dffpointer + 10,k,x); return update_coeff(s,x,-nn); >>; if m0 = 0 then << if symbolic !*tracefps then write "PS does not exist"; return failed>>; if m0 > 0 then << m1 := {}; foreach xx in denr do if type_rational rhs xx then m1 := append(list (rhs xx +1),m1); m1 := min m1; if m1 < 0 then << dff(dffpointer + 10) := df2 := x^(-m1)*dff(dffpointer); if symbolic !*tracefps then << write "working with ",x^(-m1),"*f"; write "=> f :=" ,df2 >>; s := constantre(sub(k=k+m1,cap_r),sub(k=k+m1,leadcoeff), dffpointer + 10,k,x); return update_coeff(s,x,m1) >>; >>; % { m1 >= 0 } s := 0; s0 := 0; for i:=0:m0 do << if i > m1 then dff(dffpointer + i) := df(dff(dffpointer + i-1),x); c0 := limit(dff(dffpointer + i),x,0); if (not numberp c0 and part(c0,0) = limit) then << write "Could not find the limit of: " ,dff(dffpointer + i),",",x,",",0; rederr("problem using limit operator") >> else << c0 := c0/factorial (i); s := s + c0*x^i ; if symbolic !*tracefps then write " S = ",s; >> >>; return (s); end; endmodule; module ratalgo; % rational algorithm for FPS package algebraic procedure complexapart(f,x); begin scalar !*factor,!*complex; on factor,complex; x := pf(f,x); off factor,complex; return x; end; algebraic procedure ratalgo(p,x); begin scalar afp,tt,s,ak,d,c,j,ss; afp := complexapart(p,x); s:= 0; ak := 0; if symbolic !*tracefps then write " Rational Algorithm applied"; foreach tt in afp do << if freeof(tt,x) then s := s + tt else << d := 1/tt; if not polynomq(d,x) then << if symbolic !*tracefps then write " Rational Algorithm failed"; s := -1 >>; if not (s = -1) then << d := d/lcof(d,x); c := d * tt; j := deg(d,x); d := expt(d,1/j); if not polynomq(d,x) then << if symbolic !*tracefps then write " Rational Algorithm failed"; afp := {}; d :=12; s := -1 >>; if d = x then s := s + c/d^j else << ss := lcof(d,x); d := d /ss; c := c / ss; xk := x -d; c:= c*(-1)^j/xk ^j; ak := ak + c*simplify_factorial(factorial(j + k -1)/ factorial(k)/factorial(j-1))/xk^k; >> >> >> >>; if s = -1 then return (-1); return s := s + infsum(ak*x^k,k,0,infinity) end; symbolic procedure fastpart(x,n); reval nth( x,n +1); flag ('(fastpart fastlength),'opfn); symbolic procedure fastlength(x); length (x) -1; symbolic procedure auxcopy(u); if pairp u then cons (auxcopy car u, auxcopy cdr u) else u; % for XR if getd 'print_format then print_format('(pochhammer u v),'(!( u !) !_ v)); endmodule; end;