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;