module patches; % Patches to correct problems in REDUCE 3.8.
% Author: Anthony C. Hearn.
% Copyright (c) 2004, 2005, 2006, 2007 Anthony C. Hearn. All Rights Reserved.
global '(patch!-date!* patch!-url!-list!*);
patch!-date!* := "11-Jan-2007";
patch!-url!-list!* :=
'("http://reduce-algebra.com/support/patches/patches.fsl");
% Bugs fixed by these patches.
% 26 Jun 04. With rounded arithmetic, solving some linear equation
% problems could lead to a catastrophic error.
% 8 Jul 04. Some non-zero integrals (e.g., int(e^(a^(1/3)*x)*sin x,x))
% returned zero.
% 5 Aug 04. Using RLFI with latex on could lead to invalid operator errors.
% 2 Sep 04. In rare circumstances, floating point conversion could give
% an extraneous error.
% 6 Sep 04. With rational on, some non-zero factorizations could produce
% a zero coefficient (e.g., on rational;
% factorize(r^((1/4*n^2 - 1/4*n + 1)/(n - 1)));).
% 28 Sep 04. Some integrals would not return a closed form solution
% with algint on that would with algint off
% (e.g., int(sqrt(x-1)/(sqrt x*(x-1)),x)).
% 10 Dec 04. With dfprin on, some products and sums printed incorrectly.
% 31 Jan 05. Some integrals involving square roots could run forever.
% 12 Feb 05. SOLVE could produce a spurious recursive loop (e.g.,
% solve((4*e^(y^3/3)*cte+2x^2+y^3+3)/e^(y^3/3),y)).
% 20 Apr 05. int(e^(-a^(1/4)*(-1)^(1/4)*x),x); terminated with an error.
% 2 May 05. int(e^(-a^(1/4)*(-1)^(1/4)*i*x)*b+(1/4)*e^(-a^(1/4)
% *(-1)^(1/4)*i*x)*x,x); terminated with an error.
% 22 May 05. Some integrals, e.g., int(e^((3sqrt 5+1)*x)*(sqrt 5+1)
% +e^((3sqrt 5-1)*x)*(sqrt 5+1),x), never completed.
% 30 May 05. SOLVE could produce a spurious "Zero divisor" error
% (e.g., solve({log tan(y/2),y+1/x},{x,y})).
% 4 Oct 05. DEG did not work with rational coefficients (e.g.,
% deg(x**3/a-x/5+1/4,x)).
% 5 Oct 05. Some SOLVE calculations could give a spurious "Zero Divisor"
% error (e.g., ex0:= sqrt(a^2-y^2); solve((-log(( - x + a + y)
% /ex0) + log((x + a + y)/ex0) + x - (a^2 - y^2)/ex0),y);
% 16 Nov 05. System errors could occur with rounded and combineexpt on.
% (E.g., on rounded,combineexpt; 0.183*e^x*t^4.39;).
% 22 Nov 05. Some definite integrals with variables other than x could
% give a wrong answer, e.g., int(e^(-y),y,0,x).
% 9 Dec 05. With combineexpt on, expressions could be dropped (e.g.,
% on combineexpt; 4*e^(-3*h/2) - 3*h*e^(-h) + 2*e^(-h)).
% 4 Feb 06. Setcrackflags() was not set in crack, but needed to be.
% 20 Feb 06. The rule for df(Jacobidn(~u,~m),~u) was wrong.
% 21 Feb 06. Evaluating some integrals could suppress the printing of
% the results.
% 22 Feb 06. Some sub evaluations could include superfluous terms like
% x = (x^(1/7))^7.
% 23 May 06. Derivatives and integrals of matrices were not computed.
% 18 Aug 06. After nospur, some traces were still evaluated.
% 29 Sep 06. With dfprint on, derivatives of integrals would print in a
% truncated form.
% 11 Jan 07. With rounded arithmetic and factor on, a non-numeric
% argument error could occur.
% Alg declarations.
fluid '(sublist!*);
patch alg;
% 16 Nov 05, 9 Dec 05.
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 null domainp u and null red u
and (z2 := kernels u) and null cdr z2
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 if domainp numr z2 then y := 1
else <<y := lcoeffgcd cdr comfac numr z2;
if not fixp y then 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;
% 22 Feb 06.
symbolic procedure subeval0 u;
begin scalar x,y,z,ns;
while cdr u do <<if not eqcar(car u,'equal) then x := car u . x
else if not(cadar u = (y := reval caddar u))
then x := {caar u,cadar u,y} . x;
u := cdr u>>;
if null x then return car u else u := nconc(reversip x,u);
if u member sublist!* then return mk!*sq !*p2q mksp('sub . u,1)
else sublist!* := u . sublist!*;
if null(u and cdr u)
then rederr "SUB requires at least 2 arguments";
(while cdr u do
<<x := reval car u;
if getrtype x eq 'list then u := append(cdr x,cdr u)
else <<if not eqexpr x then errpri2(car u,t);
y := cadr x;
if null getrtype y then y := !*a2kwoweight y;
if getrtype caddr x then ns := (y . caddr x) . ns
else z := (y . caddr x) . z;
u := cdr u>>>>) where !*evallhseqp=nil;
x := aeval car u;
return subeval1(append(ns,z),x)
end;
symbolic procedure subsubf(l,expn);
begin scalar x,y;
for each j in l do if car j neq (y := prepsq!* simp!* cdr j)
then x := (car j . y) . x;
l := reversip x;
if null l then return expn;
y := nil;
for each j in cddr expn do
if (x := assoc(j,l)) then <<y := x . y; l := delete(x,l)>>;
expn := sublis(l,car expn)
. for each j in cdr expn collect subsublis(l,j);
if null y then return expn;
expn := aconc!*(for each j in reversip!* y
collect list('equal,car j,aeval cdr j),expn);
return if l then subeval expn
else mk!*sq !*p2q mksp('sub . expn,1)
end;
% 23 May 06.
symbolic procedure reval1(u,v);
(begin scalar x,y;
if null u then return nil
else if stringp u then return u
else if fixp u
then return if flagp(dmode!*,'convert) then reval2(u,v) else u
else if atom u
then if null subfg!* then return u
else if idp u and (x := get(u,'avalue))
then if u memq varstack!* then recursiveerror u
else <<varstack!* := u . varstack!*;
return if y := get(car x,'evfn)
then apply2(y,u,v)
else reval1(cadr x,v)>>
else nil
else if not idp car u
then errpri2(u,t)
else if car u eq '!*sq
then return if caddr u and null !*resimp
then if null v then u else prepsqxx cadr u
else reval2(u,v)
else if flagp(car u,'remember) then return rmmbreval(u,v)
else if flagp(car u,'opfn) then return reval1(opfneval u,v)
else if x := get(car u,'psopfn)
then <<u := apply1(x,cdr u);
if x := get(x,'cleanupfn) then u := apply2(x,u,v);
return u>>
else if arrayp car u then return reval1(getelv u,v);
return if x := getrtype u then
if y := get(x,'evfn) then apply2(y,u,v)
else rerror(alg,101,
list("Missing evaluation for type",x))
else if not atom u
and not atom cdr u
and (y := getrtype cadr u)
and null(y eq 'list and cddr u)
and (x := get(y,'aggregatefn))
and (not(x eq 'matrixmap) or flagp(car u,'matmapfn))
and not flagp(car u,'boolean)
and not !*listargs and not flagp(car u,'listargp)
then apply2(x,u,v)
else reval2(u,v)
end) where varstack!* := varstack!*;
symbolic procedure getrtype2 u;
begin scalar x;
return if (x := get(car u,'rtype)) and (x := get(x,'rtypefn))
then apply1(x,cdr u)
else if x := get(car u,'rtypefn) then apply1(x,cdr u)
else if flagp(car u,'matmapfn) and cdr u
and getrtype cadr u eq 'matrix
then 'matrix
else nil
end;
endpatch;
patch arith;
% 2 Sep 04.
symbolic procedure read!:num(n);
if fixp n then make!:ibf(n, 0)
else if not(numberp n or stringp n) then bflerrmsg 'read!:num
else begin integer j,m,sign; scalar ch,u,v,l,appear!.,appear!/;
j := m := 0;
sign := 1;
u := v := appear!. := appear!/ := nil;
l := explode n;
loop: ch := car l;
if digit ch then << u := ch . u; j := j + 1 >>
else if ch eq '!. then << appear!. := t; j := 0 >>
else if ch eq '!/ then << appear!/ := t; v := u; u := nil >>
else if ch eq '!- then sign := -1
else if ch memq '(!E !D !B !e !d !b) then go to jump;
if l := cdr l then goto loop else goto make;
jump: while l := cdr l do
<<if digit(ch := car l) or ch eq '!-
then v := ch . v >>;
l := reverse v;
if car l eq '!- then m := - compress cdr l
else m:= compress l;
make: u := reverse u;
v := reverse v;
if appear!/ then
return conv!:r2bf(make!:ratnum(sign*compress v,compress u),
if !:bprec!: then !:bprec!: else 170);
if appear!. then j := - j else j := 0;
if sign = 1 then u := compress u else u := - compress u;
return round!:mt (decimal2internal (u, j + m), !:bprec!:)
where !:bprec!: := if !:bprec!: then !:bprec!:
else msd!: abs u
end;
endpatch;
patch crack;
setcrackflags();
endpatch;
% Defint declarations.
symbolic smacro procedure listsq(u);
for each uu in u collect simp!* uu;
patch defint;
symbolic procedure new_meijer(u);
begin scalar f,y,mellin,new_mellin,m,n,p,q,old_num,old_denom,temp,a1,
b1,a2,b2,alpha,num,denom,n1,temp1,temp2,coeff,v,var,new_var,new_y,
new_v,k;
f := prepsq simp car u;
y := caddr u;
mellin := bastab(car f,cddr f);
temp := car cddddr mellin;
var := cadr f;
if not idp VAR then RETURN error(99,'FAIL);
temp := reval algebraic(sub(x=var,temp));
mellin := {car mellin,cadr mellin,caddr mellin,cadddr mellin,temp};
temp := reduce_var(cadr u,mellin,var);
alpha := simp!* car temp;
new_mellin := cdr temp;
if car cddddr new_mellin neq car cddddr mellin then
<< k := car cddddr mellin;
y := reval algebraic(sub(var=y,k));
new_y := simp y>>
else
<< new_var := car cddddr new_mellin;
new_y := simp reval algebraic(sub(x=y,new_var))>>;
n1 := addsq(alpha,'(1 . 1));
temp1 := {'expt,y,prepsq n1};
temp2 := cadddr new_mellin;
coeff := simp!* reval algebraic(temp1*temp2);
m := caar new_mellin;
n := cadar new_mellin;
p := caddar new_mellin;
q := car cdddar new_mellin;
old_num := cadr new_mellin;
old_denom := caddr new_mellin;
for i:=1 :n do
<< if old_num = nil then a1 := append(a1,{simp!* old_num })
else << a1 := append(a1,{simp!* car old_num});
old_num := cdr old_num>>;
>>;
for j:=1 :m do
<< if old_denom = nil then b1 := append(b1,{simp!* old_denom })
else << b1 := append(b1,{simp!* car old_denom});
old_denom := cdr old_denom>>;
>>;
a2 := listsq old_num;
b2 := listsq old_denom;
if a1 = nil and a2 = nil then
num := list({negsq alpha})
else if a2 = nil then num := list(append(a1,{negsq alpha}))
else
<< num := append(a1,{negsq alpha}); num := append({num},a2)>>;
if b1 = nil and b2 = nil then
denom := list({subtrsq(negsq alpha,'(1 . 1))})
else if b2 = nil then
denom := list(b1,subtrsq(negsq alpha,'(1 . 1)))
else
<< denom := list(b1,subtrsq(negsq alpha,'(1 . 1)));
denom := append(denom,b2)>>;
v := gfmsq(num,denom,new_y);
if v = 'fail then return simp 'fail
else v := prepsq subsq(v,list(prepsq new_y . y));
if eqcar(v,'meijerg) then new_v := v else new_v := simp v;
return multsq(new_v,coeff);
end;
endpatch;
patch hephys;
symbolic procedure nospur u; <<rmsubs(); !*nospurp := t; flag(u,'nospur)>>;
endpatch;
% Int declarations.
fluid '(!*purerisch !*trdint gaussiani indexlist intvar sqrt!-places!-alist
loglist !*intflag!* listofnewsqrts listofallsqrts sqrt!-intvar
basic!-listofallsqrts basic!-listofnewsqrts !*precise dmode!*
!*exp !*gcd !*keepsqrts !*limitedfactors !*mcd !*rationalize
!*structure !*uncached kord!*);
smacro procedure argof u; cadr u;
patch int;
% 8 Jul 04, 31 Jan 05, 20 Apr 05, 2 May 05.
symbolic procedure df2q p;
begin scalar n,d,w,x,y,z;
if null p then return nil ./ 1;
d:=denr lc p;
w:=red p;
while w do
<<d := multf(d,quotf(denr lc w,gcdf(d,denr lc w)));
w := red w>>;
while p do begin
w := sqrt2top lc p;
x := multf(xl2f(lpow p,zlist,indexlist),multf(numr w,d));
if null x then return (p := red p);
y := denr w;
z := quotf(x,y);
if null z
then <<z := rationalizesq(x ./ y);
if denr z neq 1
then <<d := multf(denr z,d); n := multf(denr z,n)>>;
z := numr z>>;
n := addf(n,z);
p := red p
end;
return tidy!-powersq (n ./ d)
end;
% 8 Jul 04, 22 May 05.
symbolic procedure tidy!-powersq x;
begin scalar expts,!*precise,!*keepsqrts;
!*keepsqrts := t;
x := subs2q x;
expts := find!-expts(numr x,find!-expts(denr x,nil));
if null expts then return x;
x := subsq(x,for each v in expts collect
(car v . list('expt,cadr v,cddr v)));
x := subsq(x,for each v in expts collect
(cadr v
. list('expt,car v,list('quotient,1,cddr v))));
return x
end;
symbolic procedure find!-expts(ff,l);
begin scalar w;
if domainp ff then return l;
l := find!-expts(lc ff,find!-expts(red ff, l));
ff := mvar ff;
if eqcar(ff,'sqrt)
then ff := list('expt, cadr ff,'(quotient 1 2))
else if eqcar(ff,'expt) and eqcar(caddr ff,'quotient)
and numberp caddr caddr ff
then <<w := assoc(cadr ff,l);
if null w
then <<w := cadr ff . gensym() . 1; l := w . l >>;
rplacd(cdr w,lcm(cddr w,caddr caddr ff))>>;
return l
end;
% 28 Sep 04.
symbolic procedure look_for_quad(integrand, var, zz);
begin
if (car zz = 'sqrt and listp cadr zz and caadr zz = 'plus) or
(car zz = 'expt and listp cadr zz and caadr zz = 'plus and
listp caddr zz and car caddr zz = 'quotient
and fixp caddr caddr zz)
then <<
zz := simp cadr zz;
if (cdr zz = 1) then <<
zz := cdr coeff1(prepsq zz, var, nil);
if length zz = 2 then return begin
scalar a, b;
scalar nvar, res, ss;
a := car zz; b := cadr zz;
if (depends(a,var) or depends(b,var)) then return nil;
nvar := gensym();
if !*trint then <<
prin2 "Linear shift suggested ";
prin2 a; prin2 " "; prin2 b; terpri();
>>;
integrand := subsq(integrand,
list(var . list('quotient,
list('difference,
list('expt,nvar,2),a),
b)));
integrand := multsq(integrand,
simp list('quotient,list('times,nvar,2),
b));
if !*trint then <<
prin2 "Integrand is transformed by substitution to ";
printsq integrand;
prin2 "using substitution "; prin2 var; prin2 " -> ";
printsq simp list('quotient,
list('difference,list('expt,nvar,2),a),
b);
>>;
res := integratesq(integrand, nvar, nil, nil);
ss := list(nvar . list('sqrt,list('plus,list('times,var,b),
a)));
res := subsq(car res, ss) .
subsq(multsq(cdr res, simp list('quotient,b,
list('times,nvar,2))), ss);
return res;
end
else if length zz = 3 then return begin
scalar a, b, c;
a := car zz; b := cadr zz; c:= caddr zz;
if (depends(a,var) or depends(b,var) or depends(c,var)) then
return nil;
a := simp list('difference, a,
list('times,b,b,
list('quotient,1,list('times,4,c))));
if null numr a then return nil;
b := simp list('quotient, b, list('times, 2, c));
c := simp c;
return
if minusf numr c then <<
if minusf numr a then begin
scalar !*hyperbolic;
!*hyperbolic := t;
return
look_for_invhyp(integrand,nil,var,a,b,c)
end
else look_for_asin(integrand,var,a,b,c)>>
else <<
if minusf numr a then look_for_invhyp(integrand,t,var,a,b,c)
else look_for_invhyp(integrand,nil,var,a,b,c)
>>
end
else if length zz = 5 then return begin
scalar a, b, c, d, e, nn, dd, mm;
a := car zz; b := cadr zz; c:= caddr zz;
d := cadddr zz; e := car cddddr zz;
if not(b = 0) or not(d = 0) then return nil;
if (depends(a,var) or depends(c,var)) or depends(e,var) then
return nil;
nn := numr integrand; dd := denr integrand;
if denr(mm :=quotsq(nn ./ 1, !*kk2q var)) = 1 and
even_power(numr mm, var) and even_power(dd, var) then <<
return sqrt_substitute(numr mm, dd, var);
>>;
if denr(mm :=quotsq(dd ./ 1, !*kk2q var)) = 1 and
even_power(nn, var) and even_power(numr mm, var) then <<
return sqrt_substitute(nn, multf(dd,!*kk2f var), var);
>>;
return nil;
end;
>>>>;
return nil
end;
% 21 Feb 06.
symbolic procedure simpint u;
if atom u or null cdr u or cddr u and (null cdddr u or cddddr u)
then rerror(int,1,"Improper number of arguments to INT")
else if cddr u then simpdint u
else begin scalar ans,dmod,expression,variable,loglist,oldvarstack,
!*intflag!*,!*purerisch,cflag,intvar,listofnewsqrts,
listofallsqrts,sqrtfn,sqrt!-intvar,sqrt!-places!-alist,
basic!-listofallsqrts,basic!-listofnewsqrts,coefft,
varchange,w,!*precise;
!*intflag!* := t;
variable := !*a2k cadr u;
if not(idp variable or pairp variable and numlistp cdr variable)
then <<varchange := variable . intern gensym();
if !*trint
then printc {"Integration kernel", variable,
"replaced by simple variable", cdr varchange};
variable := cdr varchange>>;
intvar := variable;
w := cddr u;
if w then rerror(int,3,"Too many arguments to INT");
listofnewsqrts:= list mvar gaussiani;
listofallsqrts:= list (argof mvar gaussiani . gaussiani);
sqrtfn := get('sqrt,'simpfn);
put('sqrt,'simpfn,'proper!-simpsqrt);
if dmode!* then
<<
if (cflag:=get(dmode!*, 'cmpxfn)) then onoff('complex, nil);
if (dmod := get(dmode!*,'dname)) then
onoff(dmod,nil)>> where !*msg := nil;
begin scalar dmode!*,!*exp,!*gcd,!*keepsqrts,!*limitedfactors,!*mcd,
!*rationalize,!*structure,!*uncached,kord!*,
ans1,badbit,denexp,erfg,nexp,oneterm;
!*keepsqrts := !*limitedfactors := t;
!*exp := !*gcd := !*mcd := !*structure := !*uncached := t;
dmode!* := nil;
if !*algint
then <<
sqrt!-intvar:=!*q2f simpsqrti variable;
if (red sqrt!-intvar) or (lc sqrt!-intvar neq 1)
or (ldeg sqrt!-intvar neq 1)
then interr "Sqrt(x) not properly formed"
else sqrt!-intvar:=mvar sqrt!-intvar;
basic!-listofallsqrts:=listofallsqrts;
basic!-listofnewsqrts:=listofnewsqrts;
sqrtsave(basic!-listofallsqrts,basic!-listofnewsqrts,
list(variable . variable))>>;
coefft := (1 ./ 1);
expression := int!-simp car u;
if varchange
then <<depend1(car varchange,cdr varchange,t);
expression := int!-subsq(expression,{varchange})>>;
denexp := 1 ./ denr expression;
nexp := numr expression;
while not atom nexp and null cdr nexp and
not depends(mvar nexp,variable) do
<<coefft := multsq(coefft,(((caar nexp) . 1) . nil) ./ 1);
nexp := lc nexp>>;
ans1 := nil;
while nexp do begin
scalar x,zv,tmp;
if atom nexp then <<x := !*f2q nexp; nexp := nil>>
else <<x := !*t2q car nexp; nexp := cdr nexp>>;
x := multsq(x,denexp);
zv := zvars(getvariables x,zv,variable,t);
tmp := ans1;
while tmp do
<<if zv=caar tmp
then <<rplacd(car tmp,addsq(cdar tmp,x));
tmp := nil; zv := nil>>
else tmp := cdr tmp>>;
if zv then ans1 := (zv . x) . ans1
end;
if length ans1 = 1 then oneterm := t;
nexp := ans1;
ans := nil ./ 1;
badbit:=nil ./ 1;
while nexp do
<<u := cdar nexp;
if !*trdint
then <<princ "Integrate"; printsq u;
princ "with Zvars "; print caar nexp>>;
erfg := erfg!*;
ans1 := errorset!*(list('integratesq,mkquote u,
mkquote variable,mkquote loglist,
mkquote caar nexp),
!*backtrace);
erfg!* := erfg;
nexp := cdr nexp;
if errorp ans1 then badbit := addsq(badbit,u)
else <<ans := addsq(caar ans1, ans);
badbit:=addsq(cdar ans1,badbit)>>>>;
if !*trdint
then <<prin2 "Partial answer="; printsq ans;
prin2 "To do="; printsq badbit>>;
if badbit neq '(nil . 1)
then <<setkorder nil;
badbit := reordsq badbit;
ans := reordsq ans;
coefft := reordsq coefft;
if !*trdint then <<princ "Retrying..."; printsq badbit>>;
if oneterm and ans = '(nil . 1) then ans1 := nil
else ans1 := errorset!*(list('integratesq,mkquote badbit,
mkquote variable,mkquote loglist,nil),
!*backtrace);
if null ans1 or errorp ans1
then ans := addsq(ans,simpint1(badbit . variable . w))
else <<ans := addsq(ans,caar ans1);
if not smemq(variable, ans) then ans := nil ./ 1;
if cdar ans1 neq '(nil . 1)
then ans := addsq(ans,
simpint1(cdar ans1 . variable . w))
>>>>;
end;
ans := multsq(coefft,ans);
if !*trdint then << printc "Resimp and all that"; printsq ans >>;
put('int,'simpfn,'simpiden);
put('sqrt,'simpfn,sqrtfn);
<< if dmod then onoff(dmod,t);
if cflag then onoff('complex,t)>> where !*msg := nil;
oldvarstack := varstack!*;
varstack!* := nil;
ans := errorset!*(list('int!-resub,mkquote ans,mkquote
varchange),t);
put('int,'simpfn,'simpint);
varstack!* := oldvarstack;
return if errorp ans then error1() else car ans
end;
endpatch;
patch mathpr;
% 10 Dec 04, 29 Sep 06.
symbolic procedure dflayout u;
(begin
scalar op, args, w;
w := car (u := cdr u);
u := cdr u;
if smemq('int,w) then !*noarg := nil;
if !*noarg and (atom w or not get(car w, 'op)) then <<
if atom w then <<
op := w;
args := assoc(op, depl!*);
if args then args := cdr args >>
else <<
op := car w;
args := cdr w >>;
remember!-args(op, args);
w := op >>;
maprin w;
if u then <<
u := layout!-formula('!!dfsub!! . u, 0, nil);
if null u then return 'failed;
w := 1 + cddr u;
putpline((update!-pline(0, -w, caar u) . cdar u) .
((cadr u - w) . (cddr u - w))) >>
end) where !*noarg = !*noarg;
endpatch;
patch matrix;
% 26 Jun 04.
symbolic procedure sparse_backsub(exlis,varlis);
begin scalar d,z,c;
if null exlis then return nil;
d := lc car exlis;
foreach x in exlis do
begin scalar s,p,v,r;
p := lc x;
v := mvar x;
x := red x;
while not domainp x and mvar x member varlis do
<<if (c := atsoc(mvar x,z)) then
s := addf(multf(lc x,cdr c),s)
else r := addf(!*t2f lt x,r);
x := red x>>;
s := negf quotff(addf(multf(addf(r,x),d),s),p);
z := (v . s) . z;
end;
for each p in z do cdr p := cancel(cdr p ./ d);
return z
end;
symbolic procedure quotff(u,v);
if null u then nil
else (if x then x
else (if denr y = 1 then numr y
else rederr "Invalid division in backsub")
where y=rationalizesq(u ./ v))
where x=quotf(u,v);
% 23 May 06.
symbolic procedure matsm1 u;
begin scalar x,y,z; integer n;
a: if null u then return z
else if eqcar(car u,'!*div) then go to d
else if atom car u then go to er
else if caar u eq 'mat then go to c1
else if flagp(caar u,'matmapfn) and cdar u
and getrtype cadar u eq 'matrix
then x := matsm matrixmap(car u,nil)
else <<x := lispapply(caar u,cdar u);
if eqcar(x,'mat) then x := matsm x>>;
b: z := if null z then x
else if null cdr z and null cdar z then multsm(caar z,x)
else multm(x,z);
c: u := cdr u;
go to a;
c1: if not lchk cdar u then rerror(matrix,3,"Matrix mismatch");
x := for each j in cdar u collect
for each k in j collect xsimp k;
go to b;
d: y := matsm cadar u;
if (n := length car y) neq length y
then rerror(matrix,4,"Non square matrix")
else if (z and n neq length z)
then rerror(matrix,5,"Matrix mismatch")
else if cddar u then go to h
else if null cdr y and null cdar y then go to e;
x := subfg!*;
subfg!* := nil;
if null z then z := apply1(get('mat,'inversefn),y)
else if null(x := get('mat,'lnrsolvefn))
then z := multm(apply1(get('mat,'inversefn),y),z)
else z := apply2(get('mat,'lnrsolvefn),y,z);
subfg!* := x;
z := for each j in z collect for each k in j collect
<<!*sub2 := t; subs2 k>>;
go to c;
e: if null caaar y then rerror(matrix,6,"Zero divisor");
y := revpr caar y;
z := if null z then list list y else multsm(y,z);
go to c;
h: if null z then z := generateident n;
go to c;
er: rerror(matrix,7,list("Matrix",car u,"not set"))
end;
symbolic procedure matrixmap(u,v);
if flagp(car u,'matmapfn)
then matsm!*1 for each j in matsm cadr u collect
for each k in j collect simp!*(car u . mk!*sq k . cddr u)
else if flagp(car u,'matfn) then reval2(u,v)
else typerr(car u,"matrix operator");
put('matrix,'aggregatefn,'matrixmap);
flag('(int df taylor),'matmapfn);
flag('(det trace),'matfn);
endpatch;
patch poly;
% 6 Sep 04.
symbolic procedure fctrf u;
(begin scalar !*ezgcd,!*gcd,denom,x,y;
if domainp u then return list u
else if ncmp!* and not noncomfp u then ncmp!* := nil;
!*gcd := t;
if null !*limitedfactors and null dmode!* then !*ezgcd := t;
if null !*mcd
then rerror(poly,15,"Factorization invalid with MCD off")
else if null !*exp
then <<!*exp := t; u := !*q2f resimp !*f2q u>>;
if dmode!* eq '!:rn!:
then <<dmode!* := nil; alglist!* := nil . nil;
x := simp prepf u;
if atom denr x then <<denom := denr x; u := numr x>>
else denom := 1>>;
if null ncmp!*
then <<x := sf2ss u;
if homogp x
then <<if !*trfac
then prin2t
"This polynomial is homogeneous - variables scaled";
y := caaar x . listsum caaadr x;
x := fctrf1 ss2sf(car(x)
. (reverse subs0 cadr x . 1));
x := rconst(y,x);
return car x . sort!-factors cdr x>>>>;
u := fctrf1 u;
if denom
then <<alglist!* := nil . nil;
dmode!* := '!:rn!:; car u := quotf!*(car u,denom)>>;
return car u . sort!-factors cdr u
end) where !*exp = !*exp, ncmp!* = ncmp!*;
% 4 Oct 05.
symbolic procedure deg(u,kern);
<<u := simp!* u; tstpolyarg2(u,kern); numrdeg(numr u,kern)>>
where dmode!* = gdmode!*;
symbolic procedure tstpolyarg2(u,kern);
<<for each j in kernels numr u do
if j=kern then nil
else if depends(j,kern) then typerr(prepsq u,"polynomial");
for each j in kernels denr u do
if depends(j,kern) then typerr(prepsq u,"polynomial")>>;
% 11 Jan 07.
symbolic procedure rnfactor!: u;
begin scalar x,y,dmode!*; integer m,n;
x := subf(u,nil);
if not domainp denr x then return {1,(u . 1)};
y := factorf numr x;
n := car y;
dmode!* := '!:rn!:;
y := for each j in cdr y collect
<<n := n*(m := (lnc ckrn car j)**cdr j);
quotfd(car j,m) . cdr j>>;
return int!-equiv!-chk mkrn(n,denr x) . y
end;
endpatch;
patch rlfi;
put('tex,'simpfn,'simpcar);
endpatch;
% Solve declarations.
fluid '(!*cramer bareiss!-step!-size!*);
global '(assumptions);
patch solve;
% 26 Jun 04.
symbolic procedure solvelnrsys(exlis,varlis);
begin scalar w,x;
if w := solvesparsecheck(exlis,varlis) then exlis := w
else exlis := exlis . varlis;
if null !*cramer
and null errorp(x :=
errorset2{'solvebareiss,mkquote car exlis, mkquote cdr exlis}
where bareiss!-step!-size!* = if w then 4 else 2)
then exlis := car x
else exlis := solvecramer(car exlis,cdr exlis);
return solvesyspost(exlis,varlis)
end;
% 12 Feb 05, 5 Oct 05.
symbolic procedure solvesq (ex,var,mul);
begin scalar r,x;
r:= for each w in solvesq1(ex,var,mul) join
if null cadr w
or eqcar(x := prepsq caar w,'root_of)
or numr subfx(denr ex,{caadr w . x}) then {w};
if r and not domainp denr ex then
assumptions:=append(assumptions,{prepf denr ex});
return r
end;
% 5 Oct 05.
symbolic procedure subfx(u,v);
begin scalar x;
x := errorset2 {'subf,mkquote u,mkquote v};
return if errorp x then 1 ./ 1 else car x
end;
% 12 Feb 05
symbolic procedure polypeval u;
begin scalar bool,v;
v := cadr u;
u := simpcar u;
if cdr u neq 1 then return nil else u := kernels car u;
while u and null bool do
<<if v neq car u and smember(v,car u) then bool := t;
u := cdr u>>;
return null bool
end;
put('polyp,'psopfn,'polypeval);
(algebraic <<
depend(!~p,!~x);
clearrules
{root_of(~p,~x,~tg)^~n =>
sub(x=root_of(p,x,tg),
-reduct(p,x)/coeffn(p,x,deg(p,x)))^(n-deg(p,x)+1)
when fixp n and deg(p,x)>=1 and n>=deg(p,x)};
let root_of(~p,~x,~tg)^~n =>
sub(x=root_of(p,x,tg),
-reduct(p,x)/coeffn(p,x,deg(p,x))) ^ (n-deg(p,x)+1)
when polyp(p,x) and fixp n and deg(p,x)>=1 and n>=deg(p,x);
nodepend(!~p,!~x);
>>) where dmode!*=nil,!*modular=nil,!*rounded=nil,!*complex=nil;
% 30 May 05.
symbolic procedure solvenonlnrtansolve(u,x,w);
begin scalar v,s,z,r,y;
integer ar;
ar:=!!arbint;
v:=caar u;u:=prepf numr simp cdr u;
s:=solveeval{u,'tg!-};
!!arbint:=ar;
for each q in cdr s do
<<z:=reval caddr q;
z:=reval sublis(solvenonlnrtansolve1 z,z);
!!arbint:=ar;
y:=solve0({'equal,{'tan,{'quotient,V,2}},z},x);
r:=union(y,r)>>;
y := errorset2 {'subf,mkquote w,mkquote{x . 'pi}};
if null errorp y and null numr y
then <<!!arbint:=ar; r:=union(solve0({'equal,{'cos,x},-1},x),r)>>;
return t.r end;
% 5 Oct 05.
symbolic procedure check!-solns(z,ex,var);
begin scalar x,y;
if not errorp (x :=
errorset2 {'check!-solns1,mkquote z,mkquote ex,mkquote var})
then return car x
else if ex = (y := (numr simp!* prepf ex where !*reduced=t))
or errorp (x :=
errorset2 {'check!-solns1,mkquote z,mkquote y,mkquote var})
then return 'unsolved
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
or
errorp(y := errorset2 {'check!-solns2,mkquote ex,mkquote z})
then <<z := nil; x := 'unsolved>>
else if null(y := car 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;
symbolic procedure check!-solns2(ex,z);
if smemq('root_of,z) then rederr 'check!-solns
else numr subf(ex,{caadar z . mk!*sq caaar z});
endpatch;
patch specfn;
% 20 Feb 06.
algebraic (for all u,m let df(Jacobidn(u,m),u)
= -m^2 *Jacobisn(u,m)*Jacobicn(u,m));
endpatch;
patch rlisp;
!#if (member 'psl lispsystem!*)
symbolic procedure global idlist;
fluid idlist;
symbolic procedure global1 id1;
if not get(id1,'vartype) then fluid1 id1;
!#endif
endpatch;
endmodule;
end;