module patches; % Patches to correct problems in REDUCE 3.7.
% Author: Anthony C. Hearn.
% Copyright (c) 2001 Anthony C. Hearn. All Rights Reserved.
global '(patch!-date!*);
patch!-date!* := "6-Mar-2001";
% Bugs fixed by these patches.
% 28 Jun 99. Gnuplot handling on the Macintosh was not correct.
% 7 Aug 99. The evaluation of df(tan((sqrt(1-x^2)*asin acos x
% + 2*sqrt(1-x^2)*x)/x),x) did not terminate.
% 20 Oct 99. The sequence a1:=12x^2-16x+3; a2:=3x-4; off mcd;
% on combineexpt; e^(a1/a2); gave the wrong answer.
% 8 Nov 99. factorize(2*c*s*u^3*v^5-2*c*s*u^3*v +2*c*s*u*v^5-2*c*s*u*v
% -s^2*u^4*v^4+s^2*u^4+s^2*u^2*v^6-s^2*u^2*v^4-s^2*u^2*v^2
% +s^2*u^2 +s^2*v^6-s^2*v^2+u^4*v^4-u^4*v^2 -v^4+v^2) gave
% a catastrophic error.
% 9 Nov 99. Patched procedures generated a "redefined" message.
% 16 Nov 99. Some EXCALC calculations could cause a catastrophic error.
% 18 Dec 99. Integrations could give catastrophic errors because some
% kernels were not unique.
% 31 Jan 00. The sequence weight x=1,y=1; wtlevel 10; factor x; led to
% the error that x was invalid as a kernel.
% 5 Feb 00. The sequence x := mat((1,2)); sign sqrt 42; led to a
% spurious error.
% 6 Feb 00. The sequence on complex; sqrt(i*sqrt(3)-1); gave a wrong
% result.
% 10 Feb 00. Some root evaluations could lead to an error like
% <equation> invalid as scalar.
% 18 Feb 00. A sequence like m := mat((a,b),(c,d)); det sub(a=1,m);
% would cause a type mismatch error.
% 18 Apr 00. Complaints about the pattern matching limit of 5 terms
% are resolved by the addition of a variable matchlength!*,
% whose initial value of 5 can be changed as needed.
% 22 Apr 00. The RULE mechanism left spurious expressions in various
% non-local variables.
% 28 Jul 00. A sum index within a derivative was treated as an identifier
% (e.g., sum(x^n/factorial n*sub(x=0,df(cos x,x,n)),n,0,5);
% 2 Aug 00. With complex on, some factorizations seemed to run forever
% (e.g., factorize (400y^12+400y^10*z+40y^9*z^2+100y^8*z^2
% +20y^7*z^5+120y^7*z^4+20y^7*z^3+41y^6*z^4+60y^5*z^7
% +60y^5*z^5+20y^4*z^7+6y^4*z^6+20y^4*z^5
% +2y^3*z^6+9y^2*z^8+6y*z^8+z^8))
% 29 Aug 00. The sequence load_package gentran,scope; matrix a(10,10);
% on gentranopt; gentran a(1,1) ::=: a(1,1); caused a
% segmentation violation or similar error.
% 19 Sep 00. Clearing some sqrt rules could lead to a spurious
% "not found" message.
% 20 Sep 00. The sequence load_package algint;
% int(1/sqrt((2*e^c-y)/(e^c*y)),y);
% caused a catastrophic error.
% 8 Nov 00. Some sequences did not optimize completely when the SCOPE
% command "optimize" was used.
% 20 Nov 00. The sum operator did not always preserve a noncom order
% (e.g., noncom u,v; sum(u(m)*v(1-m),m,0,1);)
% 12 Dec 00. int with four arguments did not automatically load the
% defint package.
% 13 Dec 00. Some gcd calculations could produce an endless loop. E.g.,
% in on numval,rounded; y:=x^4+x3*x^3+x2*x^2+x1*x+x0;
% on fullroots; solve(y,x);
% 9 Jan 01. SOLVE did not return results in same order as the given
% variables (e.g., solve({y=x+t^2,x=y+u^2},{x,y,u,t});
% 14 Jan 01. Some resultants (e.g. resultant(p^3-3p^2-a,3p*(p-2),p))
% caused an error.
% 19 Jan 01. Some algebraic integrals could produce a catastrophic
% error when the algint package was loaded.
% 22 Jan 01. Some algebraic integrals could produce a spurious zero
% divisor message when the algint package was loaded (e.g.,
% int((sqrt((-sqrt(a^4*x^2+4)+a^2*x)/(2*x))
% *(-sqrt(a^4*x^2+4)*a^2*x-a^4*x^2-4))/(2*(a^4*x^2+4)),x))
% 23 Jan 01. Inverses of matrices containing non-commuting objects
% could be incorrect (e.g. noncom q; 1/mat((1,0,0),
% (x/p*q 1,1,0),(x*y/(2p*(p-1))*q 1*q 1,y/(p-2)*q 1,1))).
% 2 Feb 01. Some calls of SOLVE could produce a "zero divisor" error
% error (e.g., solve(sqrt x*sqrt((4x^2*x+1)/x)-1=0,x)).
% 9 Feb 01. The patched version of combine!-logs included an undefined
% macro.
% 20 Feb 01. Even with combineexpt on, expressions like a*a^x and
% e*e^(2/(2-x)) did not simplify adequately.
% 6 Mar 01. With algint loaded, some integrals would abort before
% completion (e.g., int((x^(2/3)*sqrt(sqrt(y)*sqrt(pi) + 2pi
% *y*x)*sqrt(- sqrt(y)*sqrt(pi)+2pi*y*x))/(4pi*y*x^3 - x),x)).
% Alg declarations.
fluid '(matchlength!*);
matchlength!* := 5;
flag('(matchlength!*),'share);
patch alg;
% 20 Oct 99, 20 Feb 01.
symbolic procedure exptunwind(u,v);
begin scalar x,x1,x2,y,z,z2;
a: if null v then return u;
x := caar v;
x1 := cadr x;
x2 := caddr x;
y := cdar v;
v := cdr v;
if !*combineexpt and length u=1 and null cdr(z2 := kernels u)
then u := {(({'expt,car z2,ldeg u} . 1) . lc u)};
while (z := assocp1(x1,v)) and
(z2 := simp {'plus,{'times,x2,y},{'times,caddar z,cdr z}})
and (!*combineexpt or (fixp numr z2 and fixp denr z2))
do <<if fixp numr z2 and fixp denr z2
then <<x2 := divide(numr z2,denr z2);
if car x2>0
then <<if fixp x1 then u := multf(x1**car x2,u)
else u := multpf(mksp(x1,car x2),u);
z2 := cdr x2 ./ denr z2>>;
y := numr z2>>
else y := 1;
x2 := prepsq(quotf(numr z2,y) ./ denr z2);
v := delete(z,v)>>;
if !*combineexpt and y=1 and fixp x1 then
<<while (z := assocp2(x2,v)) and cdr z=1 and fixp cadar z do
<<x1 := cadar z * x1; v := delete(z,v)>>;
if eqcar(x2,'quotient) and fixp cadr x2 and fixp caddr x2
and cadr x2<caddr x2
then <<z := nrootn(x1**cadr x2,caddr x2);
if cdr z = 1 then u := multd(car z,u)
else if car z = 1
then u := multf(formsf(x1,x2,1),u)
else <<u := multd(car z,u);
v := (list('expt,cdr z,x2) . 1) . v>>>>
else u := multf(formsf(x1,x2,y),u)>>
else u := multf(formsf(x1,x2,y),u);
go to a
end;
% 31 Jan 00.
symbolic procedure factor1(u,v,w);
begin scalar x,y,z,r;
y := lispeval w;
for each j in u do
if (x := getrtype j) and (z := get(x,'factor1fn))
then apply2(z,u,v)
else <<while eqcar(j:=reval j,'list) and cdr j do
<<r:=append(r,cddr j); j:=cadr j>>;
x := !*a2kwoweight j;
if v then y := aconc!*(delete(x,y),x)
else if not(x member y)
then msgpri(nil,j,"not found",nil,nil)
else y := delete(x,y)>>;
set(w,y);
if r then return factor1(r,v,w)
end;
% 5 Feb 00.
algebraic (let sign(sqrt ~a) => 1 when sign a=1);
% 18 Feb 00.
symbolic procedure getrtype u;
begin scalar x,y;
return
if null u then nil
else if atom u
then if not idp u then not numberp u and getrtype1 u
else if flagp(u,'share)
then if (x := eval u) eq u then nil else getrtype x
else if (x := get(u,'avalue)) and
not(car x memq '(scalar generic))
or (x := get(u,'rtype)) and (x := list x)
then if y := get(car x,'rtypefn) then apply1(y,nil)
else car x
else nil
else if not idp car u then nil
else if (x := get(car u,'avalue)) and (x := get(car x,'rtypefn))
then apply1(x,cdr u)
else if car u eq 'sub then 'yetunknowntype
else getrtype2 u
end;
symbolic procedure let3(u,v,w,b,flgg);
begin scalar x,y1,y2,z;
x := u;
if null x then <<u := 0; return errpri1 u>>
else if numberp x then return errpri1 u;
y2 := getrtype v;
if b and idp x then <<remprop(x,'rtype); remprop(x,'avalue)>>;
if (y1 := getrtype x)
then return if z := get(y1,'typeletfn)
then lispapply(z,list(x,v,y1,b,getrtype v))
else typelet(x,v,y1,b,getrtype v)
else if y2 and not(y2 eq 'yetunknowntype)
then return if z := get(y2,'typeletfn)
then lispapply(z,list(x,v,nil,b,y2))
else typelet(x,v,nil,b,y2)
else letscalar(u,v,w,x,b,flgg)
end;
% 18 Apr 00.
symbolic procedure mcharg1(u,v,w);
if null u and null v then list nil
else begin integer m,n;
m := length u;
n := length v;
if flagp(w,'nary) and m>2
then if m<=matchlength!* and flagp(w,'symmetric)
then return mchcomb(u,v,w)
else if n=2 then <<u := cdr mkbin(w,u); m := 2>>
else return nil;
return if m neq n then nil
else if flagp(w,'symmetric) then mchsarg(u,v,w)
else if mtp v then list pair(v,u)
else mcharg2(u,v,list nil,w)
end;
% 19 Sep 00.
symbolic procedure letscalar(u,v,w,x,b,flgg);
begin
if not atom x
then if not idp car x then return errpri2(u,'hold)
else if car x eq 'df
then if null letdf(u,v,w,x,b) then nil
else return nil
else if getrtype car x
then return let2(reval x,v,w,b)
else if not get(car x,'simpfn)
then <<redmsg(car x,"operator");
mkop car x;
return let3(u,v,w,b,flgg)>>
else nil
else if null b and null w
then <<remprop(x,'avalue);
remprop(x,'rtype);
remflag(list x,'antisymmetric);
remprop(x,'infix);
remprop(x,'kvalue);
remflag(list x,'linear);
remflag(list x,'noncom);
remprop(x,'op);
remprop(x,'opmtch);
remprop(x,'simpfn);
remflag(list x,'symmetric);
wtl!* := delasc(x,wtl!*);
if flagp(x,'opfn)
then <<remflag(list x,'opfn); remd x>>;
rmsubs();
return nil>>;
if eqcar(x,'expt) and caddr x memq frlis!*
then letexprn(u,v,w,!*k2q x,b,flgg)
else if eqcar(x,'sqrt)
then let2({'expt,cadr x,'(quotient 1 2)},v,w,
if b then b else 'sqrt);
x := simp0 x;
return if not domainp numr x then letexprn(u,v,w,x,b,flgg)
else errpri1 u
end;
symbolic procedure setk1(u,v,b);
begin scalar x,y,z,!*uncached;
!*uncached := t;
if atom u
then <<if null b
then <<if not get(u,'avalue)
then msgpri(nil,u,"not found",nil,nil)
else remprop(u,'avalue);
return nil>>
else if (x:= get(u,'avalue)) then put!-avalue(u,car x,v)
else put!-avalue(u,'scalar,v);
return v>>
else if not atom car u
then rerror(alg,25,"Invalid syntax: improper assignment");
u := car u . revlis cdr u;
if null b or b eq 'sqrt
then <<z:=assoc(u,wtl!*);
if not(y := get(car u,'kvalue))
or not (x := assoc(u,y))
then <<if null z and null(b eq 'sqrt) then
msgpri(nil,u,"not found",nil,nil)>>
else put(car u,'kvalue,delete(x,y));
if z then wtl!*:=delasc(u,wtl!*);
return nil>>
else if not (y := get(car u,'kvalue))
then put!-kvalue(car u,nil,u,v)
else <<if x := assoc(u,y)
then <<updoldrules(u,v); y := delasc(car x,y)>>;
put!-kvalue(car u,y,u,v)>>;
return v
end;
% 2 Feb 01.
symbolic procedure simprad(u,n);
if !*reduced then multsq(radfa(numr u,n),invsq radfa(denr u,n))
else begin scalar iflag,x,y,z;
if !*rationalize then <<
y:=list(denr u,1);
u:=multf(numr u, exptf(denr u,n-1)) ./ 1 >>
else y := radf(denr u,n);
if n=2 and minusf numr u
then <<iflag := t; x := radf(negf numr u,n)>>
else x := radf(numr u,n);
z := simp list('quotient,retimes cdr x, retimes cdr y);
if domainp numr z and domainp denr z
then z := multsq(mkrootsq(prepf numr z,n),
invsq mkrootsq(prepf denr z,n))
else <<if iflag
then <<iflag := nil;
z := negsq z>>;
z := mkrootsq(prepsq z,n)>>;
z := multsq(multsq(if !*precise and evenp n
then car x ./ 1
else car x ./ 1, 1 ./ car y), z);
if iflag then z := multsq(z,mkrootsq(-1,2));
return z
end;
symbolic procedure radfa(u,n);
begin scalar x,y;
x := fctrf u;
if numberp car x then x := append(zfactor car x,cdr x)
else x := (car x ./ 1) . cdr x;
y := 1 ./ 1;
for each j in x do y := multsq(y,radfb(car j,cdr j,n));
return y
end;
symbolic procedure radfb(u,m,n);
begin scalar x,y;
x := radf(u,n);
y := exptf(car x,m) ./ 1;
return multsq(exptsq(mkrootlsq(cdr x,n),m),y)
end;
% 20 Feb 01.
symbolic procedure reval2(u,v);
if v or null !*combineexpt or dmode!* then !*q2a1(simp!* u,v)
else !*q2a1((simp!* u where !*mcd = nil),v);
endpatch;
% Algint declarations.
fluid '(!*noacn !*structure !*tra !*trmin gaussiani intvar sqrtflag);
fluid '(!*pvar listofallsqrts listofnewsqrts);
global '(modevalcount);
patch algint;
% 20 Sep 00.
symbolic procedure algebraiccase(expression,zlist,varlist);
begin
scalar rischpart,deriv,w,firstterm;
scalar sqrtflag,!*structure;
sqrtflag:=t;
sqrtsave(listofallsqrts,listofnewsqrts,list(intvar . intvar));
rischpart:= errorset!*(list('doalggeom,mkquote expression),
!*backtrace);
newplace list (intvar.intvar);
if atom rischpart
then <<
if !*tra then printc "Inner integration failed";
deriv:=nil ./ 1;
rischpart:=deriv >>
else
if atom car rischpart
then <<
if !*tra or !*trmin then
printc "The 'logarithmic part' is not elementary";
return (nil ./ 1) . expression >>
else <<
rischpart:=car rischpart;
deriv:=!*diffsq(rischpart,intvar) where sqrtflag=nil;
if !*tra or !*trmin then <<
printc "Inner working yields";
printsq rischpart;
printc "with derivative";
printsq deriv >> >>;
deriv:=!*addsq(expression,negsq deriv);
if null numr deriv
then return rischpart . (nil ./ 1);
if null involvesq(deriv,intvar)
then return !*addsq(rischpart,
!*multsq(deriv,((mksp(intvar,1) .* 1) .+ nil) ./ 1))
. (nil ./ 1);
varlist:=getvariables deriv;
zlist:=findzvars(varlist,list intvar,intvar,nil);
varlist:=setdiff(varlist,zlist);
firstterm:=simp!* car zlist;
w:=sqrt2top !*multsq(deriv,invsq !*diffsq(firstterm,intvar));
if null involvesq(w,intvar)
then return !*addsq(rischpart,!*multsq(w,firstterm)) . (nil ./ 1);
if !*noacn then interr "Testing only logarithmic code";
deriv:=transcendentalcase(deriv,intvar,nil,zlist,varlist);
return !*addsq(car deriv, rischpart) . cdr deriv
end;
% 22 Jan 01, 9 Feb 01.
symbolic procedure combine!-logs(coef,logarg);
begin
scalar ans,dencoef,parts,logs,lparts,!*rationalize,trueimag;
!*rationalize:=t;
coef:=simp!* coef;
if null numr logarg then return coef;
parts:=split!-real!-imag numr coef;
if null numr cdr parts then return multsq(coef,logarg);
dencoef:=multf(denr coef,denr logarg);
if !*tra then <<
printc "attempting to find 'real' form for";
mathprint list('times,list('plus,prepsq car parts,
list('times,prepsq cdr parts,'i)),
prepsq logarg) >>;
logarg:=numr logarg;
logs:= 1 ./ 1;
while pairp logarg do <<
if ldeg logarg neq 1 then interr "what a log";
if atom mvar logarg then interr "what a log";
if car mvar logarg neq 'log then interr "what a log";
logs:=!*multsq(logs,
!*exptsq(simp!* cadr mvar logarg,lc logarg));
logarg:=red logarg >>;
logs:=rationalizesq logs;
ans:=multsq(!*multsq(car parts,logs),1 ./ dencoef);
lparts:=split!-real!-imag numr logs;
if numr difff(denr cdr lparts,intvar)
then interr "unexpected denominator";
lparts:=!*multsq(denr cdr lparts ./ 1,car lparts) . cdr lparts;
if not onep denr car lparts then interr "unexpected denominator";
trueimag:=quotsq(addf(!*exptf(numr car lparts,2),
!*exptf(numr cdr lparts,2)) ./ 1,
!*exptf(denr logs,2) ./ 1);
if numr diffsq(trueimag,intvar)
then ans:=!*addsq(ans,
!*multsq(gaussiani ./ multf(2,dencoef),
!*multsq(simplogsq trueimag,cdr parts)));
trueimag:=!*multsq(car lparts,!*invsq(numr cdr lparts ./ 1));
if numr diffsq(trueimag,intvar)
then ans:=!*addsq(ans,!*multsq(!*multsq(cdr parts,1 ./ dencoef),
!*k2q list('atan,prepsq!* trueimag)));
return ans;
end;
% 6 Mar 01.
symbolic procedure modevalvar v;
begin scalar w;
if atom v
then <<if (w := get(v,'modvalue)) then return w;
put(v,'modvalue,modevalcount);
modevalcount := modevalcount+1;
return modevalcount-1>>
else if car v neq 'sqrt
then <<if !*tra then <<princ "Unexpected algebraic:"; print v>>;
error1()>>
else if numberp cadr v then return (mksp(v,1) .* 1) .+ nil;
w := modeval(!*q2f simp cadr v,!*pvar);
w := assoc(w,listofallsqrts);
if w then return cdr w else return 'failed
end;
endpatch;
% Excalc declarations.
global '(basisforml!* detm!* indxl!* metricd!* metricu!*);
smacro procedure ldpf u;
caar u;
smacro procedure lowerind u;
list('minus,u);
patch excalc;
% 16 Nov 99.
symbolic procedure mkmetric u;
begin scalar x,y,z,okord;
putform(list(cadr u,nil,nil),0);
put(cadr u,'indxsymmetries,
'((lambda (indl) (tot!-sym!-indp
(evlis '((nth indl 1)
(nth indl 2)))))));
put(cadr u,'indxsymmetrize,
'((lambda (indl) (symmetrize!-inds '(1 2) indl))));
flag(list cadr u,'covariant);
okord := kord!*;
kord!* := basisforml!*;
x := simp!* caddr u;
y := indxl!*;
metricu!* := t;
for each j in indxl!* do
<<for each k in y do
setk(list(cadr u,lowerind j,lowerind k),0);
y := cdr y>>;
for each j on partitsq(x,'basep) do
if ldeg ldpf j = 2 then
setk(list(cadr u,lowerind cadr mvar ldpf j,
lowerind cadr mvar ldpf j),
mk!*sq lc j)
else
setk(list(cadr u,lowerind cadr mvar ldpf j,
lowerind cadr mvar lc ldpf j),
mk!*sq multsq(lc j,1 ./ 2));
kord!* := okord;
x := for each j in indxl!* collect
for each k in indxl!* collect
simpindexvar list(cadr u,lowerind j,lowerind k);
z := subfg!*;
subfg!* := nil;
y := lnrsolve(x,generateident length indxl!*);
subfg!* := z;
metricd!* := mkasmetric x;
metricu!* := mkasmetric y;
detm!* := mk!*sq detq x
end;
endpatch;
% Ezgcd declarations.
fluid '(image!-set reduced!-degree!-lclst unlucky!-case);
symbolic smacro procedure polyzerop u; null u;
patch ezgcd;
% 8 Nov 99.
symbolic procedure ezgcdf(u,v);
begin scalar kordx,x;
kordx := kord!*;
x := errorset2{'ezgcdf1,mkquote u,mkquote v};
if null errorp x then return first x;
setkorder kordx;
return gcdf1(u,v)
end;
symbolic procedure poly!-gcd(u,v);
begin scalar !*exp,z;
if polyzerop u then return poly!-abs v
else if polyzerop v then return poly!-abs u
else if u=1 or v=1 then return 1;
!*exp := t;
if quotf1(u,v) then z := v
else if quotf1(v,u) then z := u
else if !*gcd then z := gcdlist list(u,v)
else z := 1;
return poly!-abs z
end;
symbolic procedure gcdlist3(l,onestep,vlist);
begin
scalar unlucky!-case,image!-set,gg,gcont,l1,w,w1,w2,
reduced!-degree!-lclst,p1,p2;
l1:=for each p in l collect p . ezgcd!-comfac p;
l:=for each c in l1 collect
quotfail1(car c,comfac!-to!-poly cdr c,
"Content divison in GCDLIST3 failed");
gcont:=gcdlist for each c in l1 collect cddr c;
if domainp gcont then if not(gcont=1)
then errorf "GCONT has numeric part";
l := sort(for each p in l collect poly!-abs p,function ordp);
w := nil;
while l do <<
w := car l . w;
repeat l := cdr l until null l or not(car w = car l)>>;
l := reversip w;
w := nil;
if null cdr l then return multf(gcont,car l);
if domainp (gg:=car (l:=sort(l,function degree!-order))) then
return gcont;
if ldeg gg=1 then <<
if division!-test(gg,l) then return multf(poly!-abs gg,gcont)
else return gcont >>;
if onestep then <<
p1 := poly!-abs car l; p2 := poly!-abs cadr l;
if p1=p2 then <<
if division!-test(p1,cddr l) then return multf(p1,gcont) >>
else <<
gg := poly!-gcd(lc p1,lc p2);
w1 := multf(red p1, quotfail1(lc p2, gg,
"Division failure when just one pseudoremainder step needed"));
w2 := multf(red p2,negf quotfail1(lc p1, gg,
"Division failure when just one pseudoremainder step needed"));
w := ldeg p1 - ldeg p2;
if w > 0 then w2 := multf(w2, (mksp(mvar p2, w) .* 1) .+ nil)
else if w < 0
then w1 := multf(w1, (mksp(mvar p1, -w) .* 1) .+ nil);
gg := ezgcd!-pp addf(w1, w2);
if division!-test(gg,l) then return multf(gg,gcont) >>>>;
return gcdlist31(l,vlist,gcont,gg,l1)
end;
endpatch;
% Int declarations.
fluid '(tanlist);
patch int;
% 18 Dec 99.
symbolic procedure findtrialdivs zl;
begin scalar dlists1,args1;
for each z in zl do
if exportan z
then <<if car z eq 'tan
then <<args1 := (mksp(z,2) .* 1) .+ 1;
tanlist := (args1 ./ 1) . tanlist>>
else args1 := !*kk2f z;
dlists1 := (z . args1) . dlists1>>;
return dlists1
end;
% 12 Dec 00.
symbolic procedure simpdint u;
begin scalar low,upp,fn,var,x,y;
if length u neq 4
then rerror(int,2,"Improper number of arguments to INT");
load!-package 'defint;
fn := car u;
var := cadr u;
low := caddr u;
upp := cadddr u;
low := reval low;
upp := reval upp;
if low = upp then return nil ./ 1
else if null getd 'new_defint then nil
else if upp = 'infinity
then if low = 0
then if not smemql('(infinity unknown),
x := defint!* {fn,var})
then return simp!* x else nil
else if low = '(minus infinity)
then return mkinfint(fn,var)
else if freeof(var,low)
then if not smemql('(infinity unknown),
x := defint!* {fn,var})
and not smemql('(infinity unknown),
y := indefint!* {fn,var,low})
then return simp!* {'difference,x,y} else nil
else nil
else if upp = '(minus infinity) or low = 'infinity
then return negsq simpdint {fn,var,upp,low}
else if low = '(minus infinity)
then return
simpdint{prepsq simp{'sub,{'equal,var,{'minus,var}},fn},
var,{'minus,upp},'infinity}
else if low = 0
then if freeof(var,upp)
and not smemql('(infinity unknown),
x := indefint!* {fn,var,upp})
then return simp!* x else nil
else if freeof(var,upp) and freeof(var,low)
and not smemq('(infinity unknown),
x := indefint!* {fn,var,upp})
and not smemql('(infinity unknown),
y := indefint!* {fn,var,low})
then return simp!* {'difference,x,y};
return mkdint(fn,var,low,upp)
end;
endpatch;
patch matrix;
% 7 Aug 99.
symbolic procedure quotfexf!*1(u,v);
if null u then nil
else (if x then x
else (if denr y = 1 then numr y
else if denr (y := rationalizesq y)=1 then numr y
else rerror(matrix,11,
"Catastrophic division failure"))
where y=rationalizesq(u ./ v))
where x=quotf(u,v);
% 14 Jan 01.
algebraic procedure polyresultant(u,v,var);
begin scalar g,h,delta,beta,temp,uu,vv;
uu := coeff(u,var); vv := coeff(v,var);
if length uu<length vv
then return (-1 * polyresultant(v,u,var))
else if (notunivariatep(uu) > 0) or (notunivariatep(vv)>0)
then <<u := for i:=1:length uu sum
(if fixp part(uu,i) then part(uu,i)
else (co(0,part(uu,i))))*var^(i-1);
v := for i:=1:length vv sum
(if fixp part(vv,i) then part(vv,i)
else (co(0,part(vv,i))))*var^(i-1)>>;
g := h := 1;
while not (v=0) do
<<delta := deg(u,var) - deg(v,var);
beta := (-1)^(delta +1) * g * h^delta;
h := h*(lcof(v,var)/h)^delta;
temp := u;
u := v;
beta := co(0,1/beta);
v := pseudo_remainder(temp,v,var)*beta;
g := lcof(u,var)>>;
if deg(u,var) = 0 then u := u^delta else return 0;
let co_off; u := u; clearrules co_off;
return u
end;
% 23 Jan 01.
symbolic procedure lnrsolve(u,v);
begin scalar temp,vlhs,vrhs,ok,
!*exp,!*solvesingular;
if !*ncmp then return clnrsolve(u,v);
!*exp := t;
if asymplis!* or wtl!* then
<<temp := asymplis!* . wtl!*;
asymplis!* := wtl!* := nil>>;
vlhs := for i:=1:length car u collect intern gensym();
vrhs := for i:=1:length car v collect intern gensym();
u := car normmat augment(u,v);
v := append(vlhs,vrhs);
ok := setkorder v;
u := foreach r in u collect prsum(v,r);
v := errorset!*({function solvebareiss, mkquote u,mkquote vlhs},t);
if caar v memq {'singular,'inconsistent} then
<<setkorder ok; rerror(matrix,13,"Singular matrix")>>;
v := pair(cadr s,car s) where s = cadar v;
u := foreach j in vlhs collect
coeffrow(negf numr q,vrhs,denr q) where q = cdr atsoc(j,v);
setkorder ok;
if temp then <<asymplis!* := car temp; wtl!* := cdr temp>>;
return for each j in u collect
for each k in j collect
if temp then resimp k else cancel k;
end;
endpatch;
% Ncpoly declarations.
fluid '(!*complex !*trnc dipvars!*);
patch ncpoly;
% 9 Jan 01.
symbolic procedure nc_factsolve(s,vl,all);
begin scalar v,sb,ns,so,soa,sol,nz,w,q,z,r,abort;
v:= numr simp car vl;
ns:=for each e in s collect numr simp e;
r:=t;
while r do
<<r:=nil; s:=ns; ns:=nil;
for each e in s do if not abort then
<<e:=absf numr subf(e,sb);
while(q:=quotf(e,v)) do e:=q;
if null e then nil else
if domainp e or not(mvar e member vl) then abort:=t else
if null red e and domainp lc e then
<<w:=mvar e; sb:=(w . 0).sb; r:=t;
vl:=delete(w,vl)>>
else if not member(e,ns) then ns:=e.ns
>>;
>>;
if abort or null vl then return nil;
nc_factorize_timecheck();
if null ns and vl then
<<sol:={for each x in vl collect x.1};
goto done>>;
s:=for each e in ns collect prepf e;
if !*trnc then
<<prin2 "solving ";
prin2 length s; prin2 " polynomial equations for ";
prin2 length vl;
prin2t "variables";
for each e in s do writepri(mkquote e,'only);>>;
w:=(cdr solveeval{'list.s,'list.vl} where dipvars!*=nil);
loop:
nc_factorize_timecheck();
if null w then goto done;
so:= cdr car w; w:=cdr w; soa:=nil;
if smemq('i,so) and null !*complex then go to loop;
for each y in vl do if not smember(y,so) then
<<soa:=(y . 1) . soa; nz:=t>>;
for each y in so do
<<z:=nc_factorize_unwrap(reval caddr y,soa);
nz:=nz or z neq 0;
soa:=(cadr y . z).soa;
>>;
if not nz then goto loop;
q:=assoc(car vl,soa);
if null q or cdr q=0 then go to loop;
soa := for each j in soa collect (car j . sublis(soa,cdr j));
sol := soa . sol;
if all then go to loop;
done:
sol:=for each s in sol collect append(sb,s);
if !*trnc then
<<prin2t "solutions:";
for each w in sol do
writepri(mkquote('list.
for each s in w collect {'equal,car s,cdr s}),'only);
prin2t "-------------------------";
>>;
return sol
end;
endpatch;
% Plot declarations.
global '(!*plotpause !*plotusepipe dirchar!* opsys!* plotcleanup!*
plotcmds!* plotcommand!* plotdir!* plotdta!* plotheader!*
tempdir!*);
patch plot;
% 28 Jun 99.
symbolic procedure init_gnuplot();
<<
!*plotpause := -1;
plotcleanup!* := {};
tempdir!* := getenv 'tmp;
if null tempdir!* then tempdir!* := getenv 'temp;
dirchar!* := "/";
plotcommand!* := "gnuplot";
opsys!* := assoc('opsys, lispsystem!*);
if null opsys!* then opsys!* := 'unknown
else opsys!* := cdr opsys!*;
if getenv "gnuplot" then plotdir!* := getenv "gnuplot"
else if null plotdir!* and not (opsys!* = 'unix)
then plotdir!* := get!-lisp!-directory();
if opsys!* = 'win32 then <<
plotcommand!* := "wgnuplot";
plotheader!* := "";
dirchar!* := "\";
plotdta!* := for each n in
{"gnutmp.tm1", "gnutmp.tm2", "gnutmp.tm3", "gnutmp.tm4",
"gnutmp.tm5", "gnutmp.tm6", "gnutmp.tm7", "gnutmp.tm8"}
collect gtmpnam n;
plotcleanup!* := if null tempdir!* then {"erase gnutmp.tm*"}
else {bldmsg("erase %w\gnutmp.tm*", tempdir!*)} >>
else if opsys!* = 'msdos then <<
plotheader!* := ""; % ?? "set term vga";
dirchar!* := "\";
plotdta!* := for each n in
{"gnutmp.tm1", "gnutmp.tm2", "gnutmp.tm3", "gnutmp.tm4",
"gnutmp.tm5", "gnutmp.tm6", "gnutmp.tm7", "gnutmp.tm8"}
collect gtmpnam n;
plotcmds!*:= gtmpnam "gnutmp.tm0";
plotcleanup!* := if null tempdir!* then {"erase gnutmp.tm*"}
else {bldmsg("erase %w\gnutmp.tm*", tempdir!*)} >>
else if opsys!* = 'riscos then <<
plotheader!* := "";
dirchar!* := ".";
plotdta!* := for i:=1:10 collect tmpnam();
plotcmds!*:= tmpnam();
plotcleanup!* :=
bldmsg("remove %w", plotcmds!*) .
for each f in plotdta!* collect bldmsg("remove %w", f)
>>
else if opsys!* = 'unix then <<
plotheader!* := "set term x11";
plotdta!* := for i:=1:10 collect tmpnam();
plotcmds!*:= tmpnam();
plotcleanup!* :=
bldmsg("rm %w", plotcmds!*) .
for each f in plotdta!* collect bldmsg("rm %w", f) >>
else if opsys!* = 'finder then <<
plotcommand!* := "gnuplot";
plotcmds!*:= "::::gnuplot:reduce.plt";
plotheader!* := "";
dirchar!* := ":";
plotdta!* := for each n in
{"::::gnuplot:gnutmp.tm1", "::::gnuplot:gnutmp.tm2",
"::::gnuplot:gnutmp.tm3", "::::gnuplot:gnutmp.tm4",
"::::gnuplot:gnutmp.tm5", "::::gnuplot:gnutmp.tm6",
"::::gnuplot:gnutmp.tm7", "::::gnuplot:gnutmp.tm8"}
collect gtmpnam n;
plotcleanup!* := nil >>
else <<
rederr bldmsg("gnuplot for %w not available yet", opsys!*);
plotdta!* := for i:=1:10 collect tmpnam();
plotcmds!*:= tmpnam();
plotheader!* := "set term dumb" >>;
if 'pipes member lispsystem!* then !*plotusepipe:=t
else plotcommand!* := bldmsg("%w %w", plotcommand!*, plotcmds!*);
if plotdir!* then
plotcommand!* := bldmsg("%w%w%w",
plotdir!*, dirchar!*, plotcommand!*);
nil >>;
endpatch;
patch poly;
% 7 Aug 99.
symbolic procedure rationalizesq u;
begin scalar !*structure,!*sub2,v,x;
if x := get(dmode!*,'rationalizefn) then u := apply1(x,u);
powlis!* := '(i 2 (nil . t) -1 nil) . powlis!*;
v := subs2q u;
powlis!* := cdr powlis!*;
return if domainp denr v then v
else if (x := rationalizef denr v) neq 1
then <<v := multf(numr v,x) ./ multf(denr v,x);
if null !*algint and null !*rationalize
then v := gcdchk v;
subs2q v>>
else u
end;
% 6 Feb 00.
symbolic procedure sqfrf u;
begin integer n; scalar !*gcd,units,v,w,x,y,z,!*msg,r;
!*gcd := t;
if (r := !*rounded) then
<<on rational; u := numr resimp !*f2q u>>;
n := 1;
x := mvar u;
v := gcdf(u,diff(u,x));
u := quotf(u,v);
if flagp(dmode!*,'field) and ((y := lnc u) neq 1)
then <<u := multd(!:recip y,u); v := multd(y,v)>>;
while degr(v,x)>0 do
<<w := gcdf(v,u);
if u neq w then z := (quotf(u,w) . n) . z;
v := quotf(v,w);
u := w;
n := n + 1>>;
if r then
<<on rounded;
u := numr resimp !*f2q u;
z := for each j in z
collect numr resimp !*f2q car j . cdr j>>;
if v neq 1 and assoc(v,units) then v := 1;
if v neq 1 then if n=1 then u := multf(v,u)
else if (w := rassoc(1,z)) then rplaca(w,multf(v,car w))
else if null z and not domainp v then z := {v . 1}
else if null z and ((w := rootxf(v,n)) neq 'failed)
then u := multf(w,u)
else errach {"sqfrf failure",u,n,z};
return (u . n) . z
end;
% 2 Aug 00.
symbolic procedure sqfrp u;
begin scalar !*ezgcd, dmode!*;
if null getd 'ezgcdf1 then load_package ezgcd;
!*ezgcd := t;
return domainp gcdf!*(u,diff(u,mvar u))
end;
% 13 Dec 00.
symbolic procedure gcdk(u,v);
begin scalar lclst,var,w,x;
if u=v then return u
else if domainp u or degr(v,(var := mvar u))=0 then return 1
else if ldeg u<ldeg v then <<w := u; u := v; v := w>>;
if quotf1(u,v) then return v
else if !*heugcd and (x := heu!-gcd(u,v)) then return x
else if ldeg v=1
or getd 'modular!-multicheck and modular!-multicheck(u,v,var)
or not !*mcd
then return 1;
a: w := remk(u,v);
if null w then return v
else if degr(w,var)=0 then return 1;
lclst := addlc(v,lclst);
if x := quotf1(w,lc w) then w := x
else for each y in lclst do
if atom y and not flagp(dmode!*,'field)
or not
(domainp y and (flagp(dmode!*,'field)
or ((x := get(car y,'units))
and y member (for each z in x collect car z))))
then while (x := quotf1(w,y)) do w := x;
u := v; v := prim!-part w;
if degr(v,var)=0 then return 1 else go to a
end;
% 19 Jan 01.
symbolic procedure quarticf pol;
begin scalar !*sub2,a,a2,a0,b,dsc,p,p1,p2,q,shift,var;
var := mvar pol;
p := shift!-pol pol;
a := coeffs car p;
shift := caddr p;
if cadr a then rerror(poly,16,list(pol,"not correctly shifted"))
else if cadddr a then return list(1,pol);
a2 := cddr a;
a0 := caddr a2;
a2 := car a2;
a := car a;
q := quadraticf1(a,a2,a0);
if not(q eq 'failed)
then <<a2 := car q; q := cdr q;
a := exptsq(addsq(!*k2q mvar pol,shift),2);
b := numr subs2q quotsq(addsq(multsq(!*f2q car q,a),
!*f2q cadr q),
!*f2q cadr p);
a := numr subs2q quotsq(addsq(multsq(!*f2q caddr q,a),
!*f2q cadddr q),
!*f2q cadr p);
a := quadraticf!*(a,var);
b := quadraticf!*(b,var);
return multf(a2,multf(car a,car b))
. nconc!*(cdr a,cdr b)>>
else if null !*surds or denr shift neq 1
then return list(1,pol);
shift := numr shift;
if knowndiscrimsign eq 'negative then go to complex;
dsc := powsubsf addf(exptf(a2,2),multd(-4,multf(a,a0)));
p2 := minusf a0;
if not p2 and minusf dsc then go to complex;
p1 := not a2 or minusf a2;
if not p1 then if p2 then p1 := t else p2 := t;
p1 := if p1 then 'positive else 'negative;
p2 := if p2 then 'negative else 'positive;
a := rootxf(a,2);
if a eq 'failed then return list(1,pol);
dsc := rootxf(dsc,2);
if dsc eq 'failed then return list(1,pol);
p := invsq !*f2q addf(a,a);
q := multsq(!*f2q addf(a2,negf dsc),p);
p := multsq(!*f2q addf(a2,dsc),p);
b := multf(a,exptf(addf(!*k2f mvar pol,shift),2));
a := powsubsf addf(b,q);
b := powsubsf addf(b,p);
knowndiscrimsign := p1;
a := quadraticf!*(a,var);
knowndiscrimsign := p2;
b := quadraticf!*(b,var);
knowndiscrimsign := nil;
return multf(car a,car b) . nconc!*(cdr a,cdr b);
complex:
a := rootxf(a,2);
if a eq 'failed then return list(1,pol);
a0 := rootxf(a0,2);
if a0 eq 'failed then return list(1,pol);
a2 := powsubsf addf(multf(2,multf(a,a0)),negf a2);
a2 := rootxf(a2,2);
if a2 eq 'failed then return list(1,pol);
p := addf(!*k2f mvar pol,shift);
q := addf(multf(a,exptf(p,2)),a0);
p := multf(a2,p);
a := powsubsf addf(q,p);
b := powsubsf addf(q,negf p);
knowndiscrimsign := 'negative;
a := quadraticf!*(a,var);
b := quadraticf!*(b,var);
knowndiscrimsign := nil;
return multf(car a,car b) . nconc!*(cdr a,cdr b)
end;
endpatch;
% Rlisp declarations.
fluid '(newrules!*);
patch rlisp;
% 9 Nov 99.
symbolic procedure load!-package u;
begin scalar x,y;
if stringp u then return load!-package intern compress explode2 u
else if null idp u then rederr list(u,"is not a package name")
else if memq(u,loaded!-packages!*)
then return u
else if or(atom(x:= errorset(list('evload,list('quote,list u)),
nil,!*backtrace)),
cdr x)
then rederr
list("error in loading package",u,"or package not found");
loaded!-packages!* := u . loaded!-packages!*;
x := get(u,'package);
if x then x := cdr x;
a: if null x then go to b
else if null atom get(car x,'package) then load!-package car x
else if or(atom(y := errorset(list('evload,
list('quote,list car x)),
nil,!*backtrace)),
cdr y)
then rederr list("module",car x,"of package",u,
"cannot be loaded");
x := cdr x;
go to a;
b: if (x := get(u,'patchfn))
then begin scalar !*usermode,!*redefmsg; eval list x end
end;
% 22 April 00.
symbolic procedure begin11 x;
begin scalar mode,result,newrule!*;
if cursym!* eq 'end
then if terminalp() and null !*lisp!_hook
then progn(cursym!* := '!*semicol!*, !*nosave!* := t,
return nil)
else progn(comm1 'end, return 'end)
else if eqcar((if !*reduce4 then x else cadr x),'retry)
then if programl!* then x := programl!*
else progn(lprim "No previous expression",return nil);
if null !*reduce4 then progn(mode := car x,x := cadr x);
program!* := x;
if eofcheck() then return 'c else eof!* := 0;
add2inputbuf(x,if !*reduce4 then nil else mode);
if null atom x
and car x memq '(bye quit)
then if getd 'bye
then progn(lispeval x, !*nosave!* := t, return nil)
else progn(!*byeflag!* := t, return nil)
else if null !*reduce4 and eqcar(x,'ed)
then progn((if getd 'cedit and terminalp()
then cedit cdr x
else lprim "ED not supported"),
!*nosave!* := t, return nil)
else if !*defn
then if erfg!* then return nil
else if null flagp(key!*,'ignore)
and null eqcar(x,'quote)
then progn((if x then dfprint x else nil),
if null flagp(key!*,'eval) then return nil);
if !*output and ifl!* and !*echo and null !*lessspace
then terpri();
result := errorset!*(x,t);
if errorp result or erfg!*
then progn(programl!* := list(mode,x),return 'err2)
else if !*defn then return nil;
if null !*reduce4
then if null(mode eq 'symbolic) then x := getsetvars x else nil
else progn(result := car result,
(if null result then result := mkobject(nil,'noval)),
mode := type result,
result := value result);
add2resultbuf((if null !*reduce4 then car result else result),
mode);
if null !*output then return nil
else if null(semic!* eq '!$)
then if !*reduce4 then (begin
terpri();
if mode eq 'noval then return nil
else if !*debug then prin2t "Value:";
rapply1('print,list list(mode,result))
end)
else if mode eq 'symbolic
then if null car result and null(!*mode eq 'symbolic)
then nil
else begin
terpri();
result:=
errorset!*(list('print,mkquote car result),t)
end
else if car result
then result := errorset!*(list('assgnpri,mkquote car result,
(if x then 'list . x else nil),
mkquote 'only),
t);
if null !*reduce4
then return if errorp result then 'err3 else nil
else if null(!*mode eq 'noval)
then progn(terpri(), prin2 "of type: ", print mode);
return nil
end;
endpatch;
% Roots declarations.
fluid '(rootacc!#!# rootacc!#!# !*noeqns);
patch roots;
% 10 Feb 00.
symbolic procedure root_val x;
roots x
where rootacc!#!#=p, iniprec!#=p where p=precision 0, !*msg=nil,
!*noeqns=t;
endpatch;
% Scope declarations.
global '(kvarlst prevlst varlst!*);
patch scope;
% 29 Aug 00.
symbolic procedure maxtype type;
if atom type then type
else if pairp cdr type then cadr type else car type;
% 8 Nov 00.
symbolic procedure prepmultmat(preprefixlist);
begin scalar tlcm,var,varexp,kvl,kfound,pvl,pfound,tel,ratval,ratlst,
newvarlst,hvarlst;
hvarlst:= nil;
while not null (varlst!*) do
<<var := car varlst!*; varlst!* := cdr varlst!*;
if flagp(var,'ratexp)
then
<<tlcm:=1;
remflag(list var,'ratexp);
foreach elem in get(var,'varlst!*) do
if pairp cdr elem then tlcm := lcm2(tlcm,cddr elem);
varexp:=fnewsym();
tel:=(varexp.(if tlcm = 2
then list('sqrt,var)
else list('expt,var,
if onep cdr(tel:=simpquot list(1,tlcm)) then
car tel
else
list('quotient,car tel,cdr tel))));
if assoc(var,kvarlst)
then
<<kvl:=kfound:=nil;
while kvarlst and not(kfound) do
if caar(kvarlst) eq var
then
<< kvl:=tel.kvl; kfound:=t;
pvl:=pfound:=nil; prevlst:=reverse(prevlst);
while prevlst and not(pfound) do
if cdar(prevlst) eq var
then << pvl:=cons(caar prevlst,varexp).pvl;
pfound:=t
>>
else << pvl:=car(prevlst).pvl;
prevlst:=cdr(prevlst)
>>;
if pvl then
if prevlst then prevlst:=append(reverse prevlst,pvl)
else prevlst:=pvl
>>
else
<< kvl:=car(kvarlst).kvl; kvarlst:=cdr kvarlst>>;
if kvl then
if kvarlst then kvarlst:=append(reverse kvl,kvarlst)
else kvarlst:=reverse kvl
>>
else preprefixlist:=tel.preprefixlist;
ratlst:=newvarlst:=nil;
foreach elem in get(var,'varlst!*) do
if pairp cdr elem
then
<< ratval:=divide((tlcm * cadr elem)/(cddr elem),tlcm);
ratlst:=cons(car elem,cdr ratval).ratlst;
if car(ratval)>0
then newvarlst:=cons(car elem,car ratval).newvarlst
>>
else newvarlst:=elem.newvarlst;
if ratlst
then << put(varexp,'varlst!*,reverse ratlst);
hvarlst:=varexp.hvarlst
>>;
if newvarlst
then << put(var,'varlst!*,reverse newvarlst);
hvarlst:=var.hvarlst
>>
else remprop(var,'varlst!*)
>>
else hvarlst:=var.hvarlst
>>;
varlst!* := hvarlst;
return preprefixlist
end;
endpatch;
% Solve declarations.
fluid '(!*multiplicities vars!*);
global '(multiplicities!*);
patch solve;
% 9 Jan 01.
symbolic procedure !*solvelist2solveeqlist u;
begin scalar x,y,z;
u := for each j in u collect solveorder j;
for each j in u do
<<if caddr j=0 then rerror(solve,2,"zero multiplicity")
else if null cadr j
then x := for each k in car j collect
list('equal,!*q2a k,0)
else x := for each k in pair(cadr j,car j)
collect list('equal,car k,!*q2a cdr k);
if length vars!* > 1 then x := 'list . x else x := car x;
z := (caddr j . x) . z>>;
z := sort(z,function ordp);
x := nil;
if !*multiplicities
then <<for each k in z do for i := 1:car k do x := cdr k . x;
multiplicities!* := nil>>
else <<for each k in z do << x := cdr k . x; y := car k . y>>;
multiplicities!* := 'list . reversip y>>;
return 'list . reversip x
end;
symbolic procedure solveorder u;
begin scalar v,x,y,z;
v := vars!*;
x := cadr u;
if length x<length v then v := setdiff(v,setdiff(v,x));
if null x or x=v then return u;
y := car u;
while x do <<z := (car x . car y) . z; x := cdr x; y := cdr y>>;
x := for each j in v collect
if null (y := assoc(j,z)) then errach "SOLVE confusion"
else cdr y;
return x . v . cddr u
end;
% 2 Feb 01.
symbolic procedure check!-solns(z,ex,var);
begin scalar x;
if errorp (x :=
errorset2 {'check!-solns1,mkquote z,mkquote ex,mkquote var})
then return
check!-solns1(z,(numr simp!* prepf ex where !*reduced=t),var)
else return car x
end;
symbolic procedure check!-solns1(z,ex,var);
begin scalar x,y,fv,sx,vs;
fv := freevarl(ex,var);
for each z1 in z do
fv := union(fv,union(freevarl(numr caar z1,var),
freevarl(denr caar z1,var)));
fv := delete('i,fv);
if fv then for each v in fv do
if not flagp(v,'constant) then
vs := (v . list('quotient,1+random 999,1000)) . vs;
sx := if vs then numr subf(ex,vs) else ex;
while z do
if null cadar z then <<z := nil; x := 'unsolved>>
else if
<<y := numr subf(ex,list(caadar z . mk!*sq caaar z));
null y
or fv and null(y := numr subf(sx,list(caadar z .
mk!*sq subsq(caaar z,vs))))
or null numvalue y>>
then <<x := car z . x; z := cdr z>>
else z := cdr z;
return if null x then 'unsolved else x
end;
endpatch;
% Sum declarations.
fluid '(sum_last_attempt_rules!* !*zeilberg);
patch sum;
% 28 Jul 00.
symbolic procedure freeof!-df(u, v);
if atom u then t
else if car(u) eq 'df
then freeof!-df(cadr u, v) and not smember(v,cddr u)
else freeof!-dfl(cdr u, v);
symbolic procedure freeof!-dfl(u, v);
if null u then t else freeof!-df(car u,v) and freeof!-dfl(cdr u,v);
symbolic procedure simp!-sum u;
begin scalar y;
y := cdr u;
u := car u;
if not atom y and not freeof!-df(u, car y) then
if atom y
then return !*p2f(car fkern(list('sum,u)) .* 1) ./ 1
else return sum!-df(u, y);
u := simp!* u;
return if null numr u then u
else if atom y
then !*p2f(car fkern(list('sum,prepsq u)) .* 1) ./ 1
else if !*zeilberg then gosper!*(mk!*sq u,y)
else simp!-sum0(u,y)
end;
symbolic procedure sum!-subst(u,x,a);
if u = x then a
else if atom u then u
else sum!-subst(car u, x,a) . sum!-subst(cdr u,x,a);
symbolic procedure sum!-df(u,y);
begin scalar w,z,upper,lower,dif;
if length(y) = 3 then <<
lower := cadr y;
upper := caddr y;
dif := addsq(simp!* upper, negsq simp!* lower);
if denr dif = 1 then
if null numr dif
then return simp!* sum!-subst(u, car y, upper)
else if fixp numr dif then dif := numr dif
else dif := nil
else dif := nil;
if dif and dif <= 0 then return nil ./ 1 >>;
if null dif then <<
z := 'sum . (u . y);
let sum_last_attempt_rules!*;
w:= opmtch z;
rule!-list (list sum_last_attempt_rules!*,nil);
return if w then simp w else mksq(z,1)>>;
z := nil ./ 1;
a: if dif < 0 then return z;
z := addsq(z,simp!* sum!-subst(u, car y, list('plus,lower,dif)));
dif := dif - 1;
go to a
end;
% 20 Nov 00.
symbolic procedure termlst(u,v,klst);
begin scalar x,kern,lst;
if null u then return nil
else if null klst or domainp u
then return list multsq(v,!*f2q u);
kern := car klst;
klst := cdr klst;
x := setkorder list kern;
u := reorder u;
v := reorder(numr v) ./ reorder(denr v);
while not domainp u and mvar u eq kern do <<
lst := nconc(termlst(lc u, multsq(!*p2q lpow u, v),klst),lst);
u := red u>>;
if u then lst := nconc(termlst(u,v,klst),lst);
setkorder x;
return lst
end;
endpatch;
endmodule;
end;