Artifact ad5d5a27290c2befb55d24cc70c6ee7b5908fc86618f2e199cf3f5e6fd99893d:
- Executable file
r37/packages/specfn/simplede.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: 16869) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/specfn/simplede.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: 16869) [annotate] [blame] [check-ins using]
module simplede; fluid '(!*precise 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; end;