module patches; % Patches to correct problems in current release.
% Author: Anthony C. Hearn.
% Copyright (c) 1995, 1996, 1997 RAND. All Rights Reserved.
global '(patch!-date!*);
patch!-date!* := "23 Jul 98";
% Bugs fixed by these patches.
% 30 Aug 95. With depend statements, groesolve did not always eliminate
% the required variable.
% 10 Nov 95. Default value of evallhseqp changed to on, responding to
% user complaints that the "off" setting was not natural.
% 20 Nov 95. Expressions like int((1/log 2)^x,x) would not return a
% closed form.
% 21 Nov 95. Calling SUB with the wrong number of arguments could lead
% to a confusing error.
% 22 Nov 95. The modular factorizer code only works for small moduli,
% and gave strange errors if a large modulus was used.
% 27 Nov 95. Expressions like prod(k+i,i,0,0) contained "i" in the
% result. This could also cause a problem with ztrans.
% 1 Dec 95. The integral of csc(x) was not correct due to an error
% in a pattern transformation.
% 3 Dec 95. Empty rule lists led to an error.
% 4 Dec 95. On some systems, "yesp" would not respond to uppercase
% letters.
% 10 Dec 95. Primep was unnecessarily inefficient for some ranges of
% numbers.
% 12 Dec 95. The expression solve((x-i)^2,x); did not correctly show
% root multiplicity (along we several other cases).
% 18 Dec 95. In EXCALC, an index -0 could be simplified to 0.
% 19 Dec 95. Higher order partial derivatives would not simplify
% properly in EXCALC. E.g., @(f,x,y) := x**2+y**2;
% @(f,y,y,x); @(f,x,y,y); @(f,x,x,y); @(f,x,y,x);
% 27 Dec 95. REPART and IMPART didn't always compute the correct value.
% 28 Dec 95. On rare occasions, SOLVE could try to convert a rounded
% number to modular form.
% 5 Jan 96. Integrals involving i, e.g., int(i*e**(i*x)*cos(x),x)
% could generate a spurious "zero divisor" message.
% 6 Jan 96. With ROUNDED on, SOLVE could return erroneous solutions
% for problems involving elementary function expressions.
% 10 Jan 96. With ROUNDED on and SPECFN loaded, erf(-1.0) returned 0.
% 11 Jan 96. The sequence depend y,x; df(df(y,x)+3*df(y,x,2),x);
% did not recognize the dependency and returned 0.
% 12 Jan 96. solve(y*(x^x*x-1),x) produced the spurious solution y=0.
% 16 Jan 96. With SPECFN loaded, negative gamma function values could
% occur (e.g., with binomial(5,-5)).
% 17 Jan 96. The sequence let x^2=0; (1/x)^2; returned 1/0.
% 18 Jan 96. The GC switch was not defined in the CSL version.
% 26 Jan 96. The sequence load_package solve; depend y,x;
% root_of(df(y_,x)-x*y_+2,y_,tag); lost the derivative.
% 11 Feb 96. DCOMBINE would not handle TPS expressions correctly in
% some domains.
% 15 Feb 96. Integrals like int((cos(7*x))/e**(2*i*x),x) would not
% return a closed form.
% 17 Feb 96. Some output options would not work well in CSL-based
% window systems.
% 18 Feb 96. The sequence "count:=0; equ(count:=count+1):=peter;"
% set count to 2 rather than 1.
% 19 Feb 96. With some substitutions active, DET could give a
% catastrophic error.
% 25 Feb 96. Depend/nodepend updated to handle lists in first argument.
% 26 Feb 96. The following example gave incorrect zero solutions:
% solve(({ex*z+y*b1-x*b2,-z*b1+ex*y+x*b3,z*b2-y*b3+ex*x}
% where ex=sqrt(-b1**2-b2**2-b3**2)),{x,y,z});
% 27 Feb 96. With small modular coefficients, Groebner calculations
% could normalize coefficients incorrectly.
% 5 Mar 96. The sequence "load_package excalc; pform o(i,j)=0;
% o(1,20);" gave an error.
% 6 Mar 96. Errors could occur with "on ezgcd" in rounded or rational
% modes.
% 11 Mar 96. INT could fail with"on nolnr" (e.g. int(-exp(x^2), x)).
% 28 Mar 96. "on complex; e^repart x;" went into a loop.
% 15 Apr 96. The sequence "vector v1,v2; index i1,i2; nospur f;
% (v1.i1*v2.i2-v1.i2*v2.i1)*(g(f,i1,i2)-g(f,i2,i1));"
% had uncollected equal terms.
% 20 Apr 96. In rare cases, a solvable system of equations would not
% produce a solution.
% 21 Apr 96. In EXCALC, indexrange 1,2,3; and indexrange {i,j,k}=
% {1,2,3}; did not produce identical orders.
% 22 Apr 96. In EXCALC, exterior products of indexed forms did not
% respect korder settings.
% 23 Apr 96. VAROPT was not declared a switch in the solve package.
% 25 Apr 96. Some expressions involving non-commuting objects did not
% reduce to lowest terms.
% 27 Apr 96. Numerical factors were not always cancelled in numerators
% and denominators of non-commuting expressions.
% 28 Apr 96. With PRI off, non-commuting expressions were printed in
% the reverse order.
% 30 Apr 96. primep of a large prime could sometimes give an error
% " Overlarge modulus".
% 3 May 96. int(1/(x**2*y0+x**2+x*y0**2+3*x*y0+x+y0**2+y0),x) was
% returned in terms of complex expressions rather than real.
% 9 May 96. With non-id kernels, SOLVE did not always return a closed
% form when one existed (e.g., solve(-2*sqrt(x)*a
% + log(sqrt(y(x)**2 - 1) + y(x))=0, y(x))).
% 19 May 96. Solutions of some equations were not complete, e.g.,
% solve((2*y^(2/3)*a*c+2*y^(2/3)*a*x+3)/(2*y^(2/3)*a)=0,y);
% 20 May 96. The sequence load numeric;fktn:=100*(x1**2-x2)^2+(1-x1)^2;
% num_min(fktn,x1=-1.2,x2=1,iterations=200); gave wrong
% results.
% 21 May 96. The sequence "let x^3=0;on factor; x^3;" did not return 0.
% 22 May 96. Some non-commuting expressions had terms in wrong order
% after the patch of 25 Apr 96,
% e.g., noncom u,v; (u(-1)*b(1) + u(-1)*v(-1))*u(1)*b(1).
% 27 May 96. df((1/x)**(2/3),x); led to an infinite loop.
% 4 Jun 96. Solutions of some nonlinear systems were not returned in
% closed form even though straightforward solutions existed.
% 15 Jun 96. Making evallhseqp on by default caused problems in NCPOLY.
% 16 Jul 96. The sequence precision 16; on rounded;
% gcd((x-3*y)*(x-y)^5,(x-3*y)^5); returned 1.
% 17 Jul 96. det mat((0,0,a),(0,b,0),(c,0,0)); had the wrong sign.
% 18 Jul 96. Ordering problems still occurred in products that
% included noncommuting objects.
% 19 Jul 96. Some switch combinations did not simplify enough.
% e.g., off precise; on div; sqrt((a-b)/a);
% 21 Jul 96. With EXP off, df(atan((( - x + 1)/x)**0.5),x); had
% unsimplified terms.
% 27 Jul 96. Ratios of expressions with non-commuting terms were not
% always reduced to lowest terms.
% 4 Aug 96. If "factor" was used, some non-commuting expressions
% evaluations would result in an infinite loop.
% 5 Aug 96. The sequence "vector p,q; q.q where q=>p;" caused q to
% lose its vector property.
% 19 Aug 96. The sequence "p 2:= sqrt(1-p^2);" resulted in a recursive
% error.
% 11 Sep 96. The sequence "load_package arnum; on complex; on arnum;
% off arnum;" lead to a nonsensical error.
% 12 Sep 96. In very rare cases, the integration code could make an
% invalid datastructure.
% 3 Nov 96. The sequence on rational; lcof(1/2,x)+y; led to a spurious
% error.
% 2 Dec 96. sub(x=5,sqrt(x+8-6*sqrt(x-1))); returned -1 instead of 1.
% 21 Dec 96. The sequence "on evallhseqp; depend y,x; off mcd;
% {df(y,x)=(sqrt(x^4-x^2*y^2-x^2+y^2)+x^2*y-y)/(x^3-x)};"
% lead to a catastrophic error.
% 22 Dec 96. The sequence "on rounded; df(sin x,x,3.0);" gave an error.
% 20 Jan 97. solve({0,x^2},x); gave {}, not {x=0}.
% 29 Jan 97. The sequence "precision 16; on rounded,complex;
% solve({-x^3,1.002*x});" gave a system error.
% 30 Jan 97. The sequence "on rational; symbolic !:quotient(2,1);" gave
% a system error.
% 31 Jan 97. The sequence "on rounded; solve({1,1.002*x});" gave a
% system error.
% 20 Feb 97. setmod 11; on modular; solve(x=y,x); did not give a
% result in lowest terms.
% 21 Feb 97. pf(1/(x-1)^2/(x-2)/m^2,x) gave a wrong result.
% 22 Feb 97. The factorization routines did not recognize balanced_mod.
% 23 Feb 97. The sequence "on rounded; solve({1.002,0.002*x});" gave a
% system error.
% 2 Apr 97. Determinants of matrices involving surds could produce a
% catastrophic error.
% 9 Apr 97. Evaluation of some non-commuting expressions could lead to
% an infinite loop.
% 9 Aug 97. Some definite integrals were not computed to sufficient
% precision.
% 1 Sep 97. Some zero-valued expressions would not reduce properly.
% e.g., 2783/1162261467*59049**((-6)/17)*
% 127173474825648610542883299603**(23/68)*
% 43498407522777568690053259202734931579**(6/17)*
% 43498407522777568690053259202734931579**((-23)/68)
% - 2783/9*7625597484987**((-1)/34)
% *127173474825648610542883299603**(1/68)*
% 43498407522777568690053259202734931579**(6/17)*
% 43498407522777568690053259202734931579**((-23)/68);
% 21 Oct 97. Depend did not clear the algebraic cache.
% 29 Nov 97. A statement "noncom e;" incorrectly applied to the
% variable e.
% 24 Dec 97. The sequence "on complex; factorize(i*c^3+3c^2-c^2*q-3i*c
% +2i*c*q+q-1);" dropped a factor of i.
% 29 Dec 97. Expressions involving "roots_of" were not always reduced
% to lowest terms (e.g., s:=solve(p:=x^5-x+1,x); sub(s,p);
% sub(s,p)).
% 30 Dec 97. Floating point numbers with the "part" operator could
% cause a memory exception (e.g., 0.5*part({1,2},1);).
% 31 Dec 97. Successive calls of "random" did not always lead to a new
%` value.
% 1 Jan 98. The sequence "let (~x)^(~y) => ee^((~y)*log(~x)) when
% ((freeof(~x,e)) and not(numberp(~y)) and (freeof(~x,ee)));
% off mcd; (1/x-1)^-1;" gave a syntax error.
% 2 Jan 98. Some expressions with NOSPUR lines did not simplify
% completely (e.g., vector p,q; nospur l; v1 := g(l,a,p);
% v2 := g(l,a,q); v1*v2;)
% 9 Jan 98. A Hodge dual using EXCALC could include a spurious "sgn"
% term (e.g., coframe o 1 = d t, o 2 = d x; # 1;).
% 17 Jan 98. "On ezgcd" doesn't work with the complex domain. (e.g.,
% on ezgcd,complex; ( - i*m*n**3 + 3*i*m*n*q**2 + 3*i*n**3*q
% - i*n*q**3 + 3*m*n**2*q - m*q**3 + n**4 - 3*n**2*q**2)
% /(n**6 + 3*n**4*q**2 + 3*n**2*q**4 + q**6);
% 18 Jan 98. An expression like "x-y:=0" would print as "y-y:=0".
% 19 Jan 98. Matrix exponents with modular on were reduced with
% respect to the modulus.
% 21 Jan 98. With defint loaded, int(sin(y),y,0,pi) did not evaluate
% correctly.
% 24 Jan 98. The sequence "load_package excalc; pform u k=l; setmod 2;
% on modular; u(-k);" resulted in u(k) rather than u(-k).
% 27 Jan 98. In rare cases, determinants in the CSL version could have
% the wrong sign.
% 31 Jan 98. With rationalize on, expressions didn't always reduce as
% expected.
% 1 Feb 98. Some "where" statements didn't simplify fully, e.g.,
% x:= (e^(12i*pi/5) - e^(8i*pi/5) + 4e^(6i*pi/5) - e^(4i*pi/5)
% - 2e^(2i*pi/5) - 1)/(16e^(6i*pi/5)); y:= {e^(~a*i*pi/~(~ b))
% => e^((a - b)/b*i*pi) when numberp a and numberp b and a>b};
% x where y;
% 4 Feb 98. The sequence "on combineexpt,rational; 1/sqrt(2)*(rt^(2/3)
% *m^(1/3)*3^(2/3))/2^(5/6);" never completed.
% 18 Feb 98. Factorization of expressions using domains with multiple
% units could omit one of those units. E.g., on complex;
% factorize(a*( - cos(th)*a - i*r));
% 19 Feb 98. A mismatch of the dimension of a coframe basis and a
% signature gave a very confusing error message.
% E.g. coframe e 0 = d t, e 1 = d x with signature -1,1,1;
% 22 Mar 98. With precise on, some expressions were not returned as an
% absolute value (e.g., sqrt(e^(2z)/((e^(2z)+1)^2*a^2))).
% 27 Mar 98. Without the algebraic integrator loaded, some integrals of
% expressions involving fractional powers were wrong (e.g.,
% a:=r**(-3/2)*(1-r)**(-1/2)).
% 31 Mar 98. A "sqfrf failure" error could occur in some cases.
% 6 Apr 98. With algint on, some integrals could cause a catastrophic
% error (e.g., int(-2x/(sqrt(2x^2+1)-2x^2+1),x)).
% 12 May 98. The sequence k := 1; limit(k,k,0); gave the wrong value 1.
% It now gives a kernel error on k.
% 20 Jul 98. Some non-commuting expressions did not reduce to lowest
% terms (e.g., noncom p,q;for all x let p(x)*q(x)=q(x)*p(x)-i;
% cos(t)*p(x)^2*q(x)^3 - q(x)^3*cos(t)*p(x)^2;)
% 23 Jul 98. Non-commuting terms in an integrand were sometimes returned
% in the wrong order.
% Alg declarations.
fluid '(!*nospurp props!*);
patch alg;
% 10 Nov 95.
on evallhseqp;
% 21 Nov 95.
symbolic procedure subeval u;
begin scalar x,y,z,ns;
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;
if ns then x := subeval2(ns,x);
return subeval1(z,x)
end;
% 3 Dec 95.
symbolic procedure rule!-list(u,type);
begin scalar v,x,y,z;
a: frasc!* := nil;
if null u or u = {{}} then return (mcond!* := nil);
mcond!* := t;
v := car u;
if idp v
then if (x := get(v,'avalue)) and car x eq 'list
then <<u := append(reverse cdadr x,cdr u); go to a>>
else typerr(v,"rule list")
else if car v eq 'list
then <<u := append(cdr v,cdr u); go to a>>
else if car v eq 'equal
then lprim "Please use => instead of = in rules"
else if not(car v eq 'replaceby) then typerr(v,"rule");
y := remove!-free!-vars cadr v;
if eqcar(caddr v,'when)
then <<mcond!* := formbool(remove!-free!-vars!* caddr caddr v,
nil,'algebraic);
z := remove!-free!-vars!* cadr caddr v>>
else z := remove!-free!-vars!* caddr v;
rule!*(y,z,frasc!*,mcond!*,type);
u := cdr u;
go to a
end;
symbolic procedure validrule1 u;
if atom u then nil
else if car u eq 'list
then if null cdr u then {{}}
else for each j in cdr u collect validrule1 j
else if car u eq 'replaceby then u
else if car u eq 'equal then 'replaceby . cdr u
else nil;
% 10 Dec 95.
symbolic procedure primep n;
if not fixp n then typerr(n,"integer")
else if n<0 then primep(-n)
else if n=1 then nil
else if n<=!*last!-prime!-in!-list!* then n member !*primelist!*
else begin scalar p;
p := !*primelist!*;
while p and remainder(n,car p) neq 0 do p := cdr p;
return if p then nil
else if isqrt n<!*last!-prime!-in!-list!* then t
else (if n>largest!-small!-modulus then general!-primep n
else small!-primep n)
where current!-modulus=current!-modulus
end;
% 26 Jan 96, 25 Feb 96, 21 Oct 97.
symbolic procedure depend u;
depend0(u,t);
symbolic procedure nodepend u;
<<rmsubs(); depend0(u,nil)>>;
symbolic procedure depend0(u, bool);
<<alglist!* := nil . nil;
while u do
begin scalar v;
v := cdr u;
while v and not rlistp car v do v := cdr v;
for each y in (if rlistp car u then cdar u else {car u}) do
begin scalar x;
x := u;
while not((x := cdr x) eq v) do depend1(y,car x,bool);
if idp y
then <<y := intern compress append(explode y,'(!! !_));
x := u;
while not((x := cdr x) eq v) do
depend1(y,car x,bool)>>
end;
u := v
end>>;
% 15 Apr 96.
symbolic procedure ordpp(u,v);
begin scalar x;
if car u eq car v then return cdr u>cdr v;
x := kord!*;
u := car u;
v := car v;
a: if null x then return ordpa(u,v)
else if u eq car x then return t
else if v eq car x then return nil;
x := cdr x;
go to a
end;
symbolic procedure ordpa(u,v);
if null u then null v
else if null v then t
else if vectorp u then if vectorp v then ordpv(u,v) else atom v
else if atom u
then if atom v
then if numberp u then numberp v and not(u<v)
else if idp v then orderp(u,v)
else numberp v
else nil
else if atom v then t
else if car u=car v then ordpa(cdr u,cdr v)
else ordpa(car u,car v);
% 25 Apr 96, 18 Jul 96, 29 Nov 97.
symbolic procedure ordp(u,v);
if null u then null v
else if null v then t
else if vectorp u then if vectorp v then ordpv(u,v) else atom v
else if atom u
then if atom v
then if numberp u then numberp v and not(u<v)
else if idp v then orderp(u,v)
else numberp v
else nil
else if atom v then t
else if car u=car v then flagp(car u,'noncom) or ordpl(cdr u,cdr v)
else if flagp(car u,'noncom)
then if flagp(car v,'noncom) then ordp(car u, car v) else t
else if flagp(car v,'noncom) then nil
else ordp(car u,car v);
symbolic procedure ordpl(u,v);
if atom u then ordp(u,v)
else if atom v then t
else if car u=car v then ordpl(cdr u,cdr v)
else ordp(car u,car v);
% 30 Apr 96.
symbolic procedure set!-modulus p; set!-general!-modulus p;
symbolic procedure set!-general!-modulus p;
if not numberp p or p=0 then current!-modulus
else begin
scalar previous!-modulus;
previous!-modulus:=current!-modulus;
current!-modulus:=p;
modulus!/2 := p/2;
if p <= largest!-small!-modulus then set!-small!-modulus p;
return previous!-modulus
end;
% 19 May 96.
% 2 Jun 96: added test for numerical denom in quotient test.
symbolic procedure simpexpt2(u,n,flg);
begin scalar m,n,x,y;
if u=1 then return 1 ./ 1;
m:=numr n;
if pairp u then <<
if car u eq 'expt
then <<n:=multsq(m:=simp caddr u,n);
if !*precise
and numberp numr m and evenp numr m
then u := list('abs,cadr u)
else u := cadr u;
return simpexpt1(u,n,flg)>>
else if car u eq 'sqrt and not !*keepsqrts
then return simpexpt2(cadr u, multsq(1 ./ 2,n),flg)
else if car u eq 'times and not !*precise
then <<x := 1 ./ 1;
for each z in cdr u do x := multsq(simpexpt1(z,n,flg),x);
return x>>
else if car u eq 'times and (y:=split!-sign cdr u) and car y
then <<x := simpexpt1(retimes append(cadr y,cddr y),n,flg);
for each z in car y do x := multsq(simpexpt1(z,n,flg),x);
return x>>
else if car u eq 'quotient
and (not !*precise or numberp caddr u and caddr u>0)
then <<if not flg and !*mcd then
return simpexpt1(prepsq simp!* u,n,t);
n := prepsq n;
return quotsq(simpexpt{cadr u,n},simpexpt{caddr u,n})>>>>;
if null flg
then <<if null(dmode!* and idp u and get(u,dmode!*))
then u := prepsq simp!* u;
return simpexpt1(u,n,t)>>
else if numberp u and zerop u then return nil ./ 1
else if not numberp m then m := prepf m;
n := prepf denr n;
if m memq frlis!* and n=1 then return list ((u . m) . 1) . 1;
if !*mcd or not numberp m or n neq 1
or atom u or denr simp!* u neq 1 then return simpx1(u,m,n)
else return mksq(u,m)
end;
% 27 May 96.
symbolic procedure mkrootsq(u,n);
if u=1 then !*d2q 1
else if n=2 and (u= -1 or u= '(minus 1)) then simp 'i
else if eqcar(u,'expt) and fixp caddr u
then exptsq(mkrootsq(cadr u,n),caddr u)
else begin scalar x,y;
if fixp u and not minusp u
and (u<factorbound!* or !*ifactor)
and (length(x := zfactor u)>1 or cdar x>1)
then return mkrootsql(x,n);
x := if n=2 then mksqrt u
else list('expt,u,list('quotient,1,n));
if y := opmtch x then return simp y
else return mksq(x,1)
end;
% 19 Jul 96.
symbolic procedure simpexpt1(u,n,flg);
begin scalar !*allfac,!*div,m,x,y;
if onep u then return 1 ./ 1;
!*allfac := t;
m := numr n;
if m=1 and denr n=1 then return simp u;
if u eq 'e and domainp denr n and not domainp m and ldeg m=1
and null red m and eqcar(mvar m,'log) then return
simpexpt1(prepsq!* simp!* cadr mvar m,lc m ./ denr n,nil);
if not domainp m or not domainp denr n
then return simpexpt11(u,n,flg);
x := simp u;
if null m
then return if null numr x then rerror(alg,14,"0**0 formed")
else 1 ./ 1;
return if null numr x
then if domainp m and minusf m
then rerror(alg,15,"Zero divisor")
else nil ./ 1
else if atom m and denr n=1 and domainp numr x
and denr x=1
then if atom numr x and m>0 then !*d2q(numr x**m)
else <<x := !:expt(numr x,m) ./ 1;
if !*mcd then resimp x else x>>
else if y := domainvalchk('expt,list(x,n)) then y
else if atom m and denr n=1
then <<if not(m<0) then exptsq(x,m)
else if !*mcd then invsq exptsq(x,-m)
else multf(expf(numr x,m),mksfpf(denr x,-m))
./ 1>>
else simpexpt11(if flg then u else prepsq!* subs2!* x,n,t)
end;
% 5 Aug 96.
symbolic procedure evalletsub2(u,v);
begin scalar !*resimp,newrule!*,oldrules!*,props!*,x,y,z,z1;
x:=car u; u:=cadr u;
for each j in x do
if eqcar(j,'replaceby) then z := j . z
else if null v and eqcar(j,'equal)
then <<lprim "Please use => instead of = in rules";
z := ('replaceby . cdr j) . z>>
else if (y := validrule j) or (y := validrule reval j)
then (y:=reverse car y) and <<rule!-list(y,t); z1 := y . z1>>
else typerr(j,"rule list");
rule!-list(z,t);
!*resimp := t;
u := errorset!*(u,nil);
!*resimp := nil;
for each j in z . z1 do rule!-list(j,nil);
for each j in oldrules!*
do if atom cdar j
then if idp cdar j
then if cdar j eq 'scalar
then let3(caar j,cadr j,nil,t,nil)
else typelet(caar j,cadr j,nil,t,cdar j)
else nil
else rule!*(car j,cadr j,caddr j,cadddr j,'old);
restore_props();
return u
end;
symbolic procedure rule!*(u,v,frasc,mcond,type);
begin scalar x;
frasc!* := frasc;
mcond!* := mcond eq t or subla(frasc,mcond);
if type and type neq 'old
then <<newrule!* := list(u,v,frasc,mcond);
if idp u
then <<if x := get(u,'rtype)
then <<props!*:= (u . ('rtype . x)) . props!*;
remprop(u,'rtype)>>;
if x := get(u,'avalue)
then <<updoldrules(x,nil);
remprop(u,'avalue)>>>>;
if v=0 and eqcar(u,'expt) and idp cadr u
and numberp caddr u
and (x := assoc(cadr u,asymplis!*))
then updoldrules(x,nil)>>;
return rule(u,v,frasc,if type eq 'old then t else type)
end;
symbolic procedure restore_props;
for each j in props!* do
if pairp cdr j then put(car j,cadr j,cddr j)
else flag({car j},cdr j);
% 19 Aug 96.
symbolic procedure put!-kvalue(u,v,w,x);
if (if eqcar(x,'!*sq) then sq_member(w,cadr x) else smember(w,x))
then recursiveerror w
else put(u,'kvalue,aconc(v,{w,x}));
% 2 Dec 96.
symbolic procedure unnest!-sqrt!-sqrt(a,b,r); nil;
% 21 Dec 96, 22 Mar 98.
symbolic procedure radf(u,n);
begin scalar ipart,rpart,x,y,z,!*gcd,!*mcd;
if null u then return list u;
!*gcd := !*mcd := t;
ipart := 1;
z := 1;
while not domainp u do
<<y := comfac u;
if car y
then <<x := divide(pdeg car y,n);
if car x neq 0
then ipart := multf(
if evenp car x
then !*p2f(mvar u .** car x)
else if !*precise
then !*p2f mksp(mvar numr
simp list('abs,if sfp mvar u
then prepf mvar u
else mvar u),
car x)
else check!-radf!-sign(!*p2f(mvar u .** pdeg car y),
!*p2f(mvar u .** car x),
n),
ipart);
if cdr x neq 0
then rpart := mkexpt(sfchk mvar u,cdr x) . rpart>>;
x := quotf(u,comfac!-to!-poly y);
u := cdr y;
if !*reduced and minusf x
then <<x := negf x; u := negf u>>;
if flagp(dmode!*,'field) then
<<y := lnc x;
if y neq 1 then <<x := quotf(x,y); z := multd(y,z)>>>>;
if x neq 1
then <<x := radf1(sqfrf x,n);
y := car x;
if y neq 1 then
<<if !*precise and evenp n
then y := !*kk2f {'abs,prepf y};
ipart := multf(y,ipart)>>;
rpart := append(rpart,cdr x)>>>>;
if u neq 1
then <<x := radd(u,n);
ipart := multf(car x,ipart);
rpart := append(cdr x,rpart)>>;
if z neq 1
then if !*numval
and (y := domainvalchk('expt,
list(!*f2q z,!*f2q !:recip n)))
then ipart := multd(!*q2f y,ipart)
else rpart := prepf z . rpart;
return ipart . rpart
end;
% 29 Jan 97.
symbolic procedure constant_expr_listp u;
if atom u
then numberp u or flagp(u,'constant) or u eq 'i and idomainp()
else constant_exprp car u and constant_expr_listp cdr u;
% 1 Sep 97.
symbolic procedure exptchksq u;
(if u=v or denr v=1 then v
else multsq(numr v ./1,1 ./denr v))
where v = (exptchk numr u ./ exptchk denr u);
% 1 Sep 97, 4 Feb 98.
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;
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
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 := lnc 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;
symbolic procedure zfactor n; zfactor1(n,t);
symbolic procedure zfactor1(n,bool);
if n<0 then ((-1) . 1) . zfactor1(-n,bool)
else if n<4 then list(n . 1)
else begin scalar primelist,factor!-list,p,qr;
primelist := !*primelist!*;
factor!-list := nil;
while primelist do
<<p := car primelist; primelist := cdr primelist;
while cdr(qr := divide(n, p)) = 0 do
<<n:= car qr; factor!-list:= add!-factor(p,factor!-list)>>;
if n neq 1 and p*p>n
then <<primelist := nil;
factor!-list:=add!-factor(n,factor!-list);
n := 1>>>>;
if n=1 then return factor!-list;
return if null bool then (n . 1) . factor!-list
else mcfactor!*(n,factor!-list)
end;
symbolic procedure mkrootsq(u,n);
% U is a prefix expression, N an integer.
% Value is a standard quotient for U**(1/N).
if u=1 then !*d2q 1
else if n=2 and (u= -1 or u= '(minus 1)) then simp 'i
else if eqcar(u,'expt) and fixp caddr u
then exptsq(mkrootsq(cadr u,n),caddr u)
else begin scalar x,y;
if fixp u and not minusp u
and (length(x :=
zfactor1(u,u<factorbound!* or !*ifactor))>1
or cdar x>1)
then return mkrootsql(x,n);
x := if n=2 then mksqrt u
else list('expt,u,list('quotient,1,n));
if y := opmtch x then return simp y
else return mksq(x,1)
end;
% 30 Dec 97.
symbolic procedure eval!-yetunknowntypeexpr(u,v);
if atom u
then ((if w then eval!-yetunknowntypeexpr(cadr w,v)
else u)
where w = get(u,'avalue))
else if car u eq '!*sq or get(car u,'dname) or car u eq '!:dn!:
then u
else ((if x and (getrtype u eq 'yetunknowntype)
then apply1(x,cdr u)
else car u . for each j in cdr u collect
eval!-yetunknowntypeexpr(j,v))
where x = get(car u,'psopfn));
% 2 Jan 98.
symbolic procedure simp!* u;
begin scalar !*asymp!*,x;
if eqcar(u,'!*sq) and caddr u and null !*resimp
then return cadr u;
x := mul!* . !*sub2;
mul!* := nil;
u:= simp u;
if !*nospurp then mul!* := union(mul!*,'(isimpq));
for each j in mul!* do u:= apply1(j,u);
mul!* := car x;
u := subs2 u;
if !*combinelogs then u := clogsq!* u;
if dmode!* eq '!:gi!: and not !*norationalgi
then u := girationalize!: u
else if !*rationalize then u := rationalizesq u;
!*sub2 := cdr x;
if !*asymp!* and !*rationalize then u := gcdchk u;
return u
end;
% 19 Jan 98.
symbolic procedure dsimp(u,v);
if numberp u then list list u
else if atom u
then (if x and subfg!* then dsimp(cadr x,v)
else if flagp(u,'share) then dsimp(lispeval u,v)
else <<flag(list u,'used!*); list list u>>)
where x= get(u,'avalue)
else if car u eq 'plus
then for each j in cdr u join dsimp(j,v)
else if car u eq 'difference
then nconc!*(dsimp(cadr u,v),
dsimp('minus . cddr u,v))
else if car u eq 'minus
then dsimptimes(list(-1,carx(cdr u,'dsimp)),v)
else if car u eq 'times then dsimptimes(cdr u,v)
else if car u eq 'quotient
then dsimptimes(list(cadr u,list('recip,carx(cddr u,'dsimp))),v)
else if not(getrtype u eq v) then list list u
else if car u eq 'recip
then list list list('!*div,carx(cdr u,'dsimp))
else if car u eq 'expt then (lambda z;
if not numberp z then errpri2(u,t)
else if z<0
then list list list('!*div,'times . nlist(cadr u,-z))
else if z=0 then list list list('!*div,cadr u,1)
else dsimptimes(nlist(cadr u,z),v))
reval_without_mod caddr u
else if flagp(car u,'noncommuting) then list list u
else if arrayp car u
then dsimp(getelv u,v)
else (lambda x; if x then dsimp(x,v)
else (lambda y; if y then dsimp(y,v)
else list list u)
opmtch revop1 u)
opmtch u;
symbolic procedure reval_without_mod u;
if dmode!* eq '!:mod!: then (reval u where dmode!* = nil)
else reval u;
% 1 Feb 98.
symbolic procedure evalletsub1(u,v);
begin scalar x,w;
x:=car u;
u:=carx(cdr u,'simpletsub);
if eqcar(x,'list) then x := cdr x else errach 'simpletsub;
w:=evalletsub2({x,{'aeval,mkquote{'aeval,u}}},v);
return if errorp w then rerror(alg,24,"Invalid simplification")
else car w
end;
flag('(aeval),'opfn);
endpatch;
patch defint;
% 21 Jan 98.
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;
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 v;
if eqcar(v,'meijerg) then new_v := v else new_v := simp v;
return multsq(new_v,coeff);
end$
symbolic procedure reduce_var(u,v,var1);
begin scalar var,m,n,coef,alpha,beta,alpha1,alpha2,expt_flag,k,temp1,
temp2,const,new_k;
var := car cddddr(v);
beta := 1;
if length var = 0 then return u . v
else
<< k := u; coef := cadddr v;
for each i in var do
<< if listp i then
<< if car i = 'expt then
<< alpha := caddr i; expt_flag := 't>>
else if car i = 'sqrt then
<< beta := 2; alpha := 1; expt_flag := 't>>
else if car i = 'times then
<< temp1 := cadr i; temp2 := caddr i;
if listp temp1 then
<< if car temp1 = 'sqrt then
<< beta := 2; alpha1 := 1; expt_flag := 't>>
else if car temp1 = 'expt and listp caddr temp1 then
<< beta := cadr cdaddr temp1;
alpha1 := car cdaddr temp1;
expt_flag := 't>>;
>>;
if listp temp2 and car temp2 = 'expt then
<< alpha2 := caddr temp2; expt_flag := 't>>;
if alpha1 neq 'nil then
alpha := reval algebraic(alpha1 + beta*alpha2)
else alpha := alpha2;
>>;
>>
else
<< if i = 'expt then << alpha := caddr var; expt_flag := 't>>;
>>;
>>;
if expt_flag = nil then return u . v
else
<< if listp alpha then << m := cadr alpha; n := caddr alpha;
n := reval algebraic(beta*n)>>
else << m := alpha; n := beta>>;
const := reval algebraic(sub(var1=1,var));
const := reval algebraic(1/(const^(n/m)));
new_k := reval algebraic(((k + 1)*n - m)/m);
coef := reval algebraic((n/m)*coef*(const)^(k+1));
var := reval algebraic(var^(n/m));
return {new_k,car v,cadr v, caddr v,coef,var}>>;
>>;
end$
endpatch;
patch dipoly;
% 27 Feb 96.
symbolic procedure bcfi a;
(if u<0 then current!-modulus+u else u)
where u=remainder(a,current!-modulus);
endpatch;
% Excalc declarations.
fluid '(fancy!-pos!* fancy!-line!*);
global '(indxl!* keepl!* detm!* basisforml!*);
global '(basisvectorl!* dimex!* metricd!* metricu!*);
symbolic macro procedure fancy!-level u;
% unwind-protect for special output functions.
{'prog,'(pos fl w),
'(setq pos fancy!-pos!*),
'(setq fl fancy!-line!*),
{'setq,'w,cadr u},
'(cond ((eq w 'failed)
(setq fancy!-line!* fl)
(setq fancy!-pos!* pos))),
'(return w)};
symbolic smacro procedure lowerind u; list('minus,u);
smacro procedure wedgeordp(u,v); worderp(u,v);
smacro procedure !*k2pf u; u .* (1 ./ 1) .+ nil;
patch excalc;
% 18 Dec 95.
symbolic procedure putform(u,v);
if atom u then <<if flagp(u,'reserved)
then <<remflag({u},'reserved);
lpri {"***Warning: reserved variable",
u,"declared exterior form"}>>;
put(u := !*a2k u,'fdegree,list !*q2f simp v);
put(u,'clearfn,'clearfdegree)>>
else begin scalar x,y; integer n;
n := length cdr u;
if (x := get(car u,'ifdegree)) and (y := assoc(n,x))
then x := delete(y,x);
put(car u,'ifdegree,if x then (n . !*q2f simp v) . x
else list(n . !*q2f simp v));
x := car u;
flag(list x,'indexvar);
put(x,'rtype,'indexed!-form);
put(x,'simpfn,'simpindexvar);
put(x,'partitfn,'partitindexvar);
put(x,'evalargfn,'revalindl);
flag(list x,'full);
put(x,'prifn,'indvarprt);
put(x,'fancy!-pprifn,'xindvarprt);
remflag('(minus),'intfn);
if null numr simp v then flag(list x,'covariant)
end;
% 19 Dec 95.
symbolic procedure partdfpow(u,v);
begin scalar x,z; integer n;
n := cdr u;
u := car u;
z := nil ./ 1;
if u eq v then z := 1 ./ 1
else if atomf u then
if x := assoc(u,keepl!*) then
begin scalar alglist!*;
z := partdfsq(simp0 cdr x,v)
end
else if ndepends(if x := get(lid u,'varlist)
then lid u . cdr x
else lid u,v)
then z := mksq(list('partdf,u,v),1)
else return nil ./ 1
else if sfp u then z := partdff(u,v)
else if car u eq '!*sq then z := partdfsq(cadr u,v)
else if x := get(car u,dfn_prop u) then
for each j in
for each k in cdr u collect partdfsq(simp k,v)
do <<if numr j then
z := addsq(multsq(j,simp
subla(pair(caar x,cdr u),cdar x)),
z);
x := cdr x>>
else if car u eq 'partdf then
if ndepends(lid cadr u,v) then
if x := partdfsplit(u,v,get('partdf,'kvalue)) then
<<z := mksq(car x,1);
for each j in cdr x do
z := partdfsq(z,j)>>
else
<<z := 'partdf . cadr u . ordn(v . cddr u);
z := if x := opmtch z then simp x
else mksq(z,1)>>
else return nil ./ 1;
if x := atsoc(u,wtl!*) then z := multpq('k!* to (-cdr x),z);
return if n=1 then z else multsq(!*t2q((u to (n-1)) .* n),z)
end;
symbolic procedure partdfsplit(u,v,k);
if null k then nil
else if cadr caar k eq cadr u and
v memq cddr caar k and
sublistp(delete(v,cddr caar k),cddr u) then
caar k . listdiff(cddr u,delete(v,cddr caar k))
else partdfsplit(u,v,cdr k);
symbolic procedure sublistp(x,y);
null x or car x member y and sublistp(cdr x,delete(car x,y));
symbolic procedure listdiff(x,y);
if null y then x
else if null x then nil
else listdiff(delete(car y,x),cdr y);
% 5 Mar 96.
symbolic procedure xindvarprt(l,p);
fancy!-level
( if not(get('expt,'infix)>p) then
fancy!-in!-brackets(
{'xindvarprt,mkquote l,0}, '!(,'!))
else
begin scalar w,x,b,s;
w:=fancy!-prefix!-operator car l;
if w eq 'failed then return w;
l := xindxlfix cdr l;
while l and w neq 'failed do
<<if b then fancy!-prin2!*("{}",0);
b := t;
if atom car l
then (if s eq '!^
then x := car l . x
else <<if s then
<<w := fancy!-print!-indexlist1(reversip x,s,nil);
fancy!-prin2!*("{}",0)>>;
x := {car l};
s := '!^>>)
else (if s eq '!_
then x := cadar l . x
else <<if s then
<<w := fancy!-print!-indexlist1(reversip x,s,nil);
fancy!-prin2!*("{}",0)>>;
x := {cadar l};
s := '!_>>);
l := cdr l>>;
w := fancy!-print!-indexlist1(reversip x,s,nil);
return w
end);
symbolic procedure xindxlfix u;
if null u then nil
else if atom car u then xindxfix car u . xindxlfix cdr u
else {'minus,xindxfix cadar u} . xindxlfix cdr u;
symbolic procedure xindxfix x;
begin scalar xx;
xx := explode x;
while xx and car xx = '!! do xx := cdr xx;
return if xx and numberp(xx := compress xx) then xx
else x;
end;
symbolic procedure !*num2id u;
if u<10 then intern cdr assoc(u,
'((0 . !0) (1 . !1) (2 . !2) (3 . !3) (4 . !4)
(5 . !5) (6 . !6) (7 . !7) (8 . !8) (9 . !9)))
else intern compress('!! . explode u);
% 21 Apr 96.
symbolic procedure indexrange u;
begin
if null eqcar(car u,'equal)
then indxl!* := mkindxl u
else for each j in u do
begin scalar names,range;
names := cadr j;
range := caddr j;
if atom names then names := list names
else if null(car names eq 'list)
then rerror(excalc,11,
"badly formed indexrangelist")
else names := cdr names;
if atom range then range := list range
else if null(car range eq 'list)
then rerror(excalc,11,
"badly formed indexrangelist")
else range := cdr range;
range := mkindxl range;
indxl!* := reverse union(range,reverse indxl!*);
for each k in names do
put(k,'indexrange,range)
end
end;
% 22 Apr 96.
symbolic procedure indexvarordp(u,v);
if not(car u eq car v) or (u memq kord!*) or (v memq kord!*) then
ordop(u,v)
else ((if boundindp(x,indxl!*) then
if boundindp(y,indxl!*) then indordlp(cdr u,cdr v)
else t
else if boundindp(y,indxl!*) then nil
else ordop(u,v))
where x = flatindxl cdr u, y = flatindxl cdr v);
% 9 Jan 98.
symbolic procedure dual0 u;
(multpfsq(mkwedge ('wedge . basisforml!*),
simpexpt list(mk!*sq(absf!* numr x ./
absf!* denr x),'(quotient 1 2))))
where x = simp!* detm!*;
% 24 Jan 98.
symbolic procedure revalind u;
begin scalar x,y,alglist!*,dmode!*;
x := subfg!*;
subfg!* := nil;
u := subst('!0,0,u);
% The above line is used to avoid the simplifaction of -0 to 0.
y := prepsq simp u;
subfg!* := x;
return y
end;
% 19 Feb 98.
symbolic procedure cofram(u,v);
begin scalar alglist!*;
rmsubs();
u := for each j in u collect
if car j eq 'equal then cdr j else list j;
putform(caar u,1);
basisforml!* := for each j in u collect !*a2k car j;
indxl!* := for each j in basisforml!* collect cadr j;
dimex!* := length u;
basisvectorl!* := nil;
if null v then
metricd!* := nlist(1,dimex!*)
else if car v eq 'signature
then if dimex!* neq length cdr v
then rerror(excalc,12,
"Dimension of coframe and metric are inconsistent.")
else metricd!* := for each j in cdr v collect aeval j;
if null v or (car v eq 'signature) then
<<detm!* := simp car metricd!*;
for each j in cdr metricd!* do
detm!* := multsq(simp j,detm!*);
detm!* := mk!*sq detm!*;
metricu!* := metricd!*:= pair(indxl!*,for each j in
pair(indxl!*,metricd!*) collect list j)>>
else mkmetric v;
if flagp('partdf,'noxpnd) then remflag('(partdf),'noxpnd);
putform('eps . indxl!*,0);
put('eps,'indxsymmetries,
list list('lambda,'(indl),list('tot!-asym!-indp,
list('evlis,mkquote for j := 1:dimex!* collect
list('nth,'indl,j)))));
put('eps,'indxsymmetrize,
list list('lambda,'(indl),list('asymmetrize!-inds,
mkquote(for j := 1: dimex!* collect j),'indl)));
flag('(eps),'covariant);
setk('eps . for each j in indxl!* collect lowerind j,1);
if null cdar u then return;
keepl!* := append(for each j in u collect
!*a2k car j . cadr j,keepl!*);
coframe1 for each j in u collect cadr j
end;
symbolic procedure wedgek2(u,v,w);
if u eq car v and null eqcar(u,'wedge)
then if fixp deg!*form u and oddp deg!*form u then nil
else multpfsq(wedgef(u . v),mksgnsq w)
else if eqcar(car v,'wedge) then wedgek2(u,cdar v,w)
else if eqcar(u,'wedge)
then multpfsq(wedgewedge(cdr u,v),mksgnsq w)
else if wedgeordp(u,car v)
then multpfsq(wedgef(u . v),mksgnsq w)
else if cdr v
then wedgepf2(!*k2pf car v,
wedgek2(u,cdr v,addf(w,multf(deg!*form u,
deg!*form car v))))
else multpfsq(wedgef list(car v,u),
mksgnsq addf(w,multf(deg!*form u,deg!*form car v)));
endpatch;
patch factor;
% 22 Feb 97.
symbolic procedure pfactor(q,p);
begin scalar user!-prime,current!-modulus,modulus!/2,r,x;
if not numberp p then typerr(p,"number")
else if not primep p then typerr(p,"prime")
else if p>largest!-small!-modulus
then rederr {p,"too large a modulus for factorization"};
user!-prime:=p;
set!-modulus p;
if domainp q or null reduce!-mod!-p lc q then
prin2t "*** Degenerate case in modular factorization";
if not (length variables!-in!-form q=1) then
return fctrfkronm q;
r:=reduce!-mod!-p q;
x := lnc r;
r :=monic!-mod!-p r;
r:=errorset!*(list('factor!-form!-mod!-p,mkquote r),t);
if not errorp r
then return x . for each j in car r
collect mod!-adjust car j . cdr j;
prin2t "****** FACTORIZATION FAILED******";
return list(1,prepf q)
end;
symbolic procedure mod!-adjust u;
% Make sure any modular numbers in u are in the right range.
if null !*balanced_mod then u else mod!-adjust1 u;
symbolic procedure mod!-adjust1 u;
if domainp u
then if fixp u then !*modular2f u
else if eqcar(u,'!:mod!:) then !*modular2f cdr u
else typerr(u,"modular number")
else lpow u .* mod!-adjust1 lc u .+ mod!-adjust1 red u;
endpatch;
% Groebner declarations.
fluid '(!*trgroeb intvdpvars!*);
patch groebner;
% 30 Aug 95.
symbolic procedure vdpvordopt (w,vars);
begin; scalar c;
vars := sort(vars,'ordop);
c := for each x in vars collect x . 0 . 0;
for each poly in w do vdpvordopt1 (poly,vars,c);
c := sort(c,function vdpvordopt2);
intvdpvars!* := for each v in c collect car v;
vars := vdpvordopt31 intvdpvars!*;
if !*trgroeb then
<<prin2 "optimized sequence of kernels: "; prin2t vars>>;
return (for each poly in w collect reorder poly) . vars;
end;
symbolic procedure vdpvordopt31 u;
begin scalar v,y;
if null u then return nil;
v := foreach x in u join
<< y := assoc(x,depl!*);
if null y or null xnp(cdr y,u) then {x} >>;
return nconc(vdpvordopt31 setdiff(u,v),v);
end;
endpatch;
patch hephys;
% 2 Jan 98.
remflag('(nospur),'flagop);
symbolic procedure nospur u; <<!*nospurp := t; flag(u,'nospur)>>;
endpatch;
% Int declarations.
fluid '(!*trint intvar lorder power!-list!* lhs!* rhs!* ulist zlist);
fluid '(!*galois !*hyperbolic !*pvar listofnewsqrts);
fluid '(!*algint !*purerisch !*trdint gaussiani sqrt!-places!-alist);
global '(!*number!* !*statistics);
patch int;
% 20 Nov 95.
symbolic procedure exportan dl;
if atom dl then nil
else begin
if eq(car dl,'tan) then return t;
nxt: if not eq(car dl,'expt) then return nil;
dl:=cadr dl;
if atom dl or not smember(intvar,dl) then return t;
go to nxt
end;
% 1 Dec 95.
symbolic procedure hfaglk(k,x);
begin
scalar kt;
if atom k or not member(x,flatten(cdr(k))) then return !*k2q k;
k := car(k) . hfaglargs(cdr(k), x);
if cadr k eq 'pi then return !*k2q k;
kt := simp list('tan, list('quotient, cadr(k), 2));
return if car(k) = 'sin
then quotqq(multsq(!*int2qm(2),kt), addsq(!*int2qm(1),
exptsq(kt,2)))
else if car(k) = 'cos
then quotqq(!*subtrq(!*int2qm(1),exptsq(kt,2)),addsq(!*int2qm(1),
exptsq(kt,2)))
else if car(k) = 'tan
then quotqq(multsq(!*int2qm(2),kt), !*subtrq(!*int2qm(1),
exptsq(kt,2)))
else if car(k) = 'cot
then quotqq(!*subtrq(!*int2qm(1),
exptsq(kt,2)),multsq(!*int2qm(2),kt))
else if car(k) = 'sec
then quotqq(addsq(!*int2qm(1), exptsq(kt,2)),
!*subtrq(!*int2qm(1),exptsq(kt,2)))
else if car(k) = 'csc
then quotqq(addsq(!*int2qm(1),exptsq(kt,2)),
multsq(!*int2qm(2),kt))
else if car(k) = 'sinh then
quotqq(!*subtrq(!*p2q mksp('expt.('e. cdr k),2),
!*int2qm(1)), multsq(!*int2qm(2),
!*p2q mksp('expt . ('e . cdr(k)),1)))
else if car(k) = 'cosh then
quotqq(addsq(!*p2q mksp('expt.('e. cdr k),2),
!*int2qm(1)), multsq(!*int2qm(2),
!*p2q mksp('expt . ('e . cdr(k)),1)))
else if car(k) = 'tanh then
quotqq(!*subtrq(!*p2q mksp('expt.('e. cdr k),2),
!*int2qm(1)), addsq(!*p2q mksp ('expt.('e.cdr(k)),2),
!*int2qm(1)))
else if car(k) = 'coth then
quotqq(addsq(!*p2q mksp('expt.('e.cdr(k)),2), !*int2qm(1)),
!*subtrq(!*p2q mksp('expt.('e . cdr k),2),!*int2qm(1)))
else if car(k) = 'acot then
!*p2q mksp(list('atan, list('quotient, 1, cadr k)),1)
else !*k2q(k);
end;
% 5 Jan 96.
symbolic procedure solve!-for!-u(rhs!*,lhs!*,ulist);
begin
top:
if null lhs!* then return ulist
else begin scalar u,lpowlhs;
lpowlhs := lpow lhs!*;
begin scalar ll,m1,chge;
ll:=maxorder(power!-list!*,zlist,0);
m1:=lorder;
while m1 do << if car ll < car m1 then
<< chge:=t; rplaca(m1,car ll) >>;
ll:=cdr ll; m1:=cdr m1 >>;
if !*trint and chge then <<
princ
"Maximum order for undetermined coefficients is reduced to ";
printc lorder >>
end;
u:=pickupu(rhs!*,lpow lhs!*,t);
if null u then
<< if !*trint then <<
printc "***** Equation for a constant to be solved:";
printsf numr lc lhs!*;
printc " = 0";
printc " ">>;
if gausselimn(numr lc lhs!*,lt lhs!*) then <<
lhs!*:=squashconstants(red lhs!*); u := t >>
else lhs!*:=red lhs!* >>
else
<< ulist:=(car u .
subs2q multsq(coefdf(lhs!*,lpowlhs),invsq cdr u)).ulist;
if !*statistics then !*number!*:=!*number!*+1;
if !*trint then <<
printc "A coefficient of numerator has been determined";
prin2 "***** U"; prin2 car u; prin2t " =";
printsq multsq(coefdf(lhs!*,lpowlhs),invsq cdr u);
printc " ">>;
lhs!*:=plusdf(lhs!*,
negdf multdfconst(cdar ulist,uterm(car u,rhs!*)))>>;
if !*trint and u
then <<printc "Terms remaining are:"; printdf lhs!*;
printc " ">>
end;
go to top
end;
% 15 Feb 96.
symbolic procedure multdfconst(x,u);
if null u or null numr x then nil
else lpow u .* subs2q multsq(x,lc u) .+ multdfconst(x,red u);
% 11 Mar 96.
symbolic procedure simpint!* u;
begin scalar x;
return if (x := opmtch('int . u)) then simp x
else simpiden('int!* . u)
end;
% 3 May 96.
symbolic procedure int!-quadterm(pol,var);
begin scalar a,b,c,discrim,kord,res,w;
kord := setkorder(var . kord!*);
c := reorder pol;
if ldeg c neq 2
then <<setkorder kord;
rerror(int,5,"Invalid polynomial in int-quadterm")>>;
a := lc c;
c := red c;
if not domainp c and mvar c = var and ldeg c = 1
then <<b := lc c; c := red c>>;
setkorder kord;
discrim := powsubsf addf(multf(b,b),multd(-4,multf(a,c)));
if null discrim then interr "discrim is zero in quadterm";
w := rootxf!*(negf discrim,2);
if not(w eq 'failed) then go to atancase;
w := rootxf!*(discrim,2);
if w eq 'failed then return list ('log . !*f2q pol);
discrim := w;
w := multpf(mksp(var,1),a);
w := addf(multd(2,w),b);
a := addf(w,discrim);
b := addf(w,negf discrim);
a := quotf(a,cdr comfac a);
b := quotf(b,cdr comfac b);
return ('log . !*f2q a) . ('log . !*f2q b) . res;
atancase:
res := ('log . !*f2q pol) . res;
a := multpf(mksp(var,1),a);
a := addf(b,multd(2,a));
a := fquotf(a,w);
return ('atan . a) . res
end;
symbolic procedure rootxf!*(u,n);
(if x eq 'failed or smemq('i,x) and not smemq('i,u)
then (rootxf(u,n) where !*surds=nil)
else x)
where x=rootxf(u,n);
% 12 Sep 96.
symbolic procedure rootextractf v;
if domainp v then v
else begin scalar u,r,c,x,p;
u := mvar v; p := ldeg v;
r := rootextractf red v;
c := rootextractf lc v;
if null c then return r
else if atom u then return (lpow v .* c) .+ r
else if car u eq 'sqrt
or car u eq 'expt and eqcar(caddr u,'quotient)
and car cdaddr u = 1 and numberp cadr cdaddr u
then <<p := divide(p,if car u eq 'sqrt then 2
else cadr cdaddr u);
if car p = 0
then return if null c then r else (lpow v .* c) .+ r
else if numberp cadr u
then <<c := multd(cadr u ** car p,c); p := cdr p>>
else <<x := simpexpt list(cadr u,car p);
if denr x = 1 then <<c := multf(numr x,c); p := cdr p>>
else p := ldeg v>>>>;
return if p=0 then addf(c,r)
else if null c then r
else ((u to p) .* c) .+ r
end;
% 27 Mar 98.
symbolic procedure look_for_quad(integrand, var, zz);
if !*algint then nil
else 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;
symbolic procedure look_for_invhyp(integrand, do_acosh, var, a, b, c);
begin
scalar newvar, res, ss, sqmn, onemth, m, n, realdom;
m := prepsq a;
n := prepsq c;
b := prepsq b;
newvar := gensym();
if do_acosh then <<
sqmn := prepsq apply1(get('sqrt, 'simpfn),
list list('quotient, n, list('minus, m)));
onemth := list('sinh, newvar);
ss := list('cosh, newvar)
>>
else <<
sqmn:= prepsq apply1(get('sqrt,'simpfn),list list('quotient,n,m));
onemth := list('cosh, newvar);
ss := list('sinh, newvar)
>>;
powlis!* := list(ss, 2, '(nil . t),
list((if do_acosh then 'plus else 'difference),
list('expt, onemth, 2),1),
nil) .
powlis!*;
integrand := subs2q
multsq(subsq(integrand,
list(var . list('difference,list('quotient,ss,sqmn),b))),
quotsq(onemth := simp onemth, simp sqmn));
if !*trint then <<
prin2 "Integrand is transformed by substitution to ";
printsq integrand;
prin2 "using substitution "; prin2 var; prin2 " -> ";
printsq simp list('difference, list('quotient, ss, sqmn), b);
>>;
realdom := not smember('(sqrt -1),integrand);
res := integratesq(integrand, newvar, nil, nil);
powlis!* := cdr powlis!*;
if !*hyperbolic then <<
ss := list(if do_acosh then 'acosh else 'asinh,
list('times,list('plus,var,b), sqmn));
>>
else <<
ss := list('times,list('plus,var,b), sqmn);
ss := if do_acosh then
subst(ss,'ss,
'(log (plus ss (sqrt (difference (times ss ss) 1)))))
else
subst(ss,'ss,'(log (plus ss (sqrt (plus (times ss ss) 1)))))
>>;
ss := list(newvar . ss);
res := sqrt2top subsq(car res, ss) .
sqrt2top subsq(quotsq(cdr res, onemth), ss);
if !*trint then <<
printc "Transforming back...";
printsq car res;
prin2 " plus a bad part of ";
printsq cdr res
>>;
if (car res = '(nil . 1)) then return nil;
if realdom and smember('(sqrt -1),res) then <<
if !*trint then print "Wrong sheet"; return nil;
>>;
return res
end;
% 6 Apr 98.
symbolic procedure mknewsqrt u;
begin scalar v,w;
if not !*galois then go to new;
v:=addf(!*p2f mksp(!*pvar,2),negf !*q2f tidysqrt simp u);
w:=errorset!*(list('afactor,mkquote v,mkquote !*pvar),t);
if atom w then go to new
else w:=car w;
if cdr w then go to notnew;
new:
w := mksqrt reval u;
listofnewsqrts:=w . listofnewsqrts;
return !*kk2f w;
notnew:
w:=car w;
v:=stt(w,!*pvar);
if car v neq 1 then errach list("Error in mknewsqrt: ",v);
w:=addf(w,multf(cdr v,(mksp(!*pvar,car v) .* -1) .+nil));
v:=sqrt2top(w ./ cdr v);
w:=quotf(numr v,denr v);
if null w
then begin scalar oldprop;
oldprop := get('sqrt,'simpfn);
put('sqrt,'simpfn,'simpsqrt);
v := simp prepsq v;
put('sqrt,'simpfn,oldprop);
if denr v = 1 then w := numr v
end;
if null w then errach list("Division failure in mknewsqrt",u);
return w
end;
% 23 Jul 98.
symbolic procedure !*multf(u,v);
begin scalar x,y;
if null u or null v then return nil
else if u=1 then return squashsqrt v
else if v=1 then return squashsqrt u
else if domainp u then return multd(u,squashsqrt v)
else if domainp v then return multd(v,squashsqrt u)
else if noncomfp u or noncomfp v then return multf(u,v);
x:=mvar u;
y:=mvar v;
if x eq y then go to c else if ordop(x,y) then go to b;
x:=!*multf(u,lc v);
y:=!*multf(u,red v);
return if null x then y
else if not domainp lc v
and mvar u eq mvar lc v
and not atom mvar u
and car mvar u memq '(expt sqrt)
then addf(!*multf(x,!*p2f lpow v),y)
else makeupsf(lpow v,x,y);
b: x:=!*multf(lc u,v);
y:=!*multf(red u,v);
return if null x then y
else if not domainp lc u
and mvar lc u eq mvar v
and not atom mvar v
and car mvar v memq '(expt sqrt)
then addf(!*multf(!*p2f lpow u,x),y)
else makeupsf(lpow u,x,y);
c: y:=addf(!*multf(list lt u,red v),!*multf(red u,v));
if eqcar(x,'sqrt)
then return addf(squashsqrt y,!*multfsqrt(x,
!*multf(lc u,lc v),ldeg u + ldeg v))
else if eqcar(x,'expt) and prefix!-rational!-numberp caddr x
then return addf(squashsqrt y,!*multfexpt(x,
!*multf(lc u,lc v),ldeg u + ldeg v));
x:=mkspm(x,ldeg u + ldeg v);
return if null x or null (u:=!*multf(lc u,lc v))
then y
else addf(multpf(x,u),y)
end;
% 23 Jul 98.
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;
!*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,denexp,badbit,nexp,oneterm,!*precise;
!*keepsqrts := !*limitedfactors := t;
!*exp := !*gcd := !*mcd := !*structure := !*uncached := t;
dmode!* := nil;
if !*algint
then <<
!*precise := t;
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 := findzvars(getvariables x,list variable,variable,nil);
begin scalar oldzlist;
while oldzlist neq zv do <<
oldzlist := zv;
foreach zz in oldzlist do
zv:=findzvars(distexp(pseudodiff(zz,variable)),
zv,variable,t)>>;
zv := sort(zv, function ordp)
end;
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>>;
ans1 := errorset!*(list('integratesq,mkquote u,
mkquote variable,mkquote loglist,
mkquote caar nexp),
!*backtrace);
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;
% 17 Feb 96.
symbolic procedure fvarpri(u,v,w);
begin integer scountr,llength,nchars; scalar explis,fvar,svar;
fortlang!* := reval fort_lang;
if not(fortlang!* memq '(fort c)) then
typerr(fortlang!*,"target language");
if not posintegerp card_no
then typerr(card_no,"FORTRAN card number");
if not posintegerp fort_width
then typerr(fort_width,"FORTRAN line width");
llength := linelength fort_width;
if stringp u
then return <<fprin2!* u;
if w eq 'only then fterpri(t);
linelength llength>>;
if eqcar(u,'!*sq) then u := prepsq!* sqhorner!* cadr u;
scountr := 0;
nchars := if fortlang!* = 'c then 999999
else ((linelength nil-spare!*)-12)*card_no;
svar := varnam!*;
fvar := if null v then (if fortlang!*='fort then svar else nil)
else car v;
if posn!*=0 and w then fortpri(fvar,u,w)
else fortpri(nil,u,w);
linelength llength
end;
% 28 Apr 96.
symbolic procedure xprinf(u,v,w);
begin
while not domainp u do <<xprint(lt u,v); u := red u; v := t>>;
if null u then return nil
else if minusf u then <<oprin 'minus; u := !:minus u>>
else if v then oprin 'plus;
if atom u then prin2!* u else maprin u
end;
symbolic procedure xprint(u,flg);
begin scalar v,w;
v := tc u;
u := tpow u;
if (w := kernlp v) and w neq 1
then <<v := quotf(v,w);
if minusf w
then <<oprin 'minus; w := !:minus w; flg := nil>>>>;
if flg then oprin 'plus;
if w and w neq 1 then <<prin2!* w; oprin 'times>>;
xprinp u;
if v neq 1
then <<oprin 'times;
if red v then prin2!* "(";
xprinf(v,nil,nil);
if red v then prin2!* ")">>
end;
symbolic procedure xprinp u;
begin
if atom car u then prin2!* car u
else if not atom caar u or caar u eq '!*sq then <<
prin2!* "(";
if not atom caar u then xprinf(car u,nil,nil)
else sqprint cadar u;
prin2!* ")" >>
else if caar u eq 'plus then maprint(car u,100)
else maprin car u;
if (u := cdr u)=1 then return nil
else if !*nat and !*eraise
then <<ycoord!* := ycoord!*+1;
if ycoord!*>ymax!* then ymax!* := ycoord!*>>
else prin2!* get('expt,'prtch);
prin2!* if numberp u and minusp u then list u else u;
if !*nat and !*eraise
then <<ycoord!* := ycoord!*-1;
if ymin!*>ycoord!* then ymin!* := ycoord!*>>
end;
% 18 Jan 98.
symbolic procedure evalvars u;
if null u then nil
else if atom car u or flagp(caar u,'intfn)
then car u . evalvars cdr u
else (caar u . revlis cdar u) . evalvars cdr u;
endpatch;
% Matrix declarations.
fluid '(bareiss!-step!-size!*);
global '(assumptions requirements);
patch matrix;
% 19 Feb 96, 27 Jan 98.
symbolic procedure bareiss!-det u;
begin scalar nu,bu,n,ok,temp,v,!*exp;
!*exp := t;
nu := matsm car u;
n := length nu;
for each x in nu do
if length x neq n then rederr "Non square matrix";
if n=1 then return caar nu;
if asymplis!* or wtl!* then
<<temp := asymplis!* . wtl!*;
asymplis!* := wtl!* := nil>>;
nu := normmat nu;
v := for i:=1:n collect intern gensym();
ok := setkorder append(v,kord!*);
car nu := foreach r in car nu collect prsum(v,r);
begin scalar powlis,powlis1;
powlis := powlis!*; powlis!* := nil;
powlis1 := powlis1!*; powlis1!* := nil;
bu := cdr sparse_bareiss(car nu,v,bareiss!-step!-size!*);
powlis!* := powlis; powlis1!* := powlis1;
end;
bu := if length bu = n then (lc car bu ./ cdr nu) else (nil ./ 1);
setkorder ok;
if temp then <<asymplis!* := car temp; wtl!* := cdr temp>>;
return resimp bu
end;
% 17 Jul 96.
symbolic procedure sparse_bareiss(u,v,k);
begin scalar p,d,w,pivs,s;
d := 1;
u := foreach f in u join if f then {!*sf2ex(f,v)};
while p := choose_pivot_rows(u,v,k,d) do
begin
u := car p; v := cadr p;
p := cddr p;
pivs := lpow car p;
u := foreach r in u join
begin
r := splitup(r,v);
r := extadd(extmult(cadr r,car p),extmult(car r,cadr p));
if null (r := subs2chkex r) then return nil;
r := innprodpex(pivs,quotexf!*(r,d));
if not evenp length pivs then r := negex r;
return {r};
end;
d := lc car p;
assumptions := 'list . mk!*sq !*f2q d .
(pairp assumptions and cdr assumptions);
p := extadd(car p,cadr p);
s := evenp length pivs;
foreach x in pivs do
w := if (s := not s) then innprodpex(delete(x,pivs),p) . w
else negex innprodpex(delete(x,pivs),p) . w;
end;
foreach f in u do
requirements := 'list . mk!*sq !*f2q !*ex2sf f .
(pairp requirements and cdr requirements);
return if u then 'inconsistent . nil
else t . foreach f in w collect !*ex2sf f;
end;
symbolic procedure choose_pivot_rows(u,v,k,d);
if null u or null v then nil else
begin scalar w,s,ss,p,x,y,rows,pivots;
w := u;
for i:=1:k do if v then v := cdr v;
while k neq 0 do
if null u then
if null v or null w or pivots then k := 0
else
<< for i:=1:k do if v then v := cdr v; s := nil; u := w>>
else if car(x := splitup(car u,v)) and
(y := if null pivots then car x
else subs2chkex extmult(car x,car pivots)) then
begin
rows := x . rows;
pivots := (if null pivots then y else quotexf!*(y,d)) . pivots;
if s then ss := not ss;
w := delete(car u,w); u := cdr u; k := k - 1;
end
else <<u := cdr u; s := not s>>;
if null pivots then return;
if remainder(length lpow car pivots,4) member {2,3} then
ss := not ss;
rows := reverse rows;
pivots := reverse pivots;
p := car rows;
foreach r in cdr rows do
p := {car(pivots := cdr pivots),
quotexf!*(extadd(extmult(cadr r,car p),
extmult(car r,cadr p)),d)};
return w . v . if ss then {negex car p,negex cadr p} else p;
end;
% 2 Apr 97.
symbolic procedure simpdet u;
begin scalar x,!*protfg;
!*protfg := t;
return if !*cramer or !*rounded or
errorp(x := errorset({'bareiss!-det,mkquote u},nil,nil))
then detq matsm carx(u,'det)
else car x
end;
symbolic procedure quotexf!*(u,v);
if null u then nil
else lpow u .* quotfexf!*1(lc u,v) .+ quotexf!*(red u,v);
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 rerror(matrix,11,"Bareiss division failure"))
where y=rationalizesq(u ./ v))
where x=quotf(u,v);
endpatch;
patch misc;
% 27 Nov 95.
symbolic procedure simp!-prod u;
begin scalar v,y,upper,lower,lower1,dif;
y := cdr u;
u := simp!* car u;
if null numr u then return (1 ./ 1)
else if atom y then return u;
if not atom cdr y then <<
lower := cadr y;
lower1 := if numberp lower then lower - 1
else list('plus,lower,-1);
upper := if not atom cddr y then caddr y else car y;
dif := addsq(simp!* upper, negsq simp!* lower);
if denr dif = 1 then
if null numr dif
then return subsq(u,list(!*a2k 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 1 ./ 1;
if atom cddr y then upper := nil>>;
v := !*a2k car y;
return simp!-prod1(u,v,y,upper,lower,lower1,dif)
end;
% 21 Feb 97.
symbolic procedure pf(u,var);
begin scalar !*exp,!*gcd,kord!*,!*limitedfactors,polypart,rfactor,
u1,u2,u3,u4,var,x,xx,y;
!*exp := !*gcd := t;
xx := updkorder var;
x := subs2 resimp simp!* u;
u1 := denr x;
if degr(u1,var) = 0 then <<setkorder xx; return list('list,u)>>;
u2 := qremsq(!*f2q numr x,!*f2q u1,var);
if caar u2 then polypart := car u2;
rfactor := 1 ./ 1;
u2 := cdr u2;
u3 := fctrf u1;
x := cdr u3;
u3 := car u3;
while not domainp u3 do
<<if mvar u3 eq var then x := (!*k2f mvar u3 . ldeg u3) . x
else <<u4 := !*p2f lpow u3;
rfactor := numr rfactor ./ multf(u4,denr rfactor);
u1 := quotf(u1,u4)>>;
u3 := lc u3>>;
if u3 neq 1 then <<rfactor := numr rfactor
./ multf(u3,denr rfactor);
u1 := quotf(u1,u3)>>;
while length x>1 do
<<u3 := exptf(caar x,cdar x);
u1 := quotf(u1,u3);
if degr(u3,var)=0
then rfactor := numr rfactor ./ multf(u3,denr rfactor)
else <<u4 := xeucl(u1,u3,var);
y := (multsq(remsq(multsq(car u4,u2),!*f2q u3,var),
rfactor) . car x)
. y;
u2 := multsq(cdr u4,u2)>>;
x := cdr x>>;
u3 := exptf(caar x,cdar x);
if u2 = (nil ./ 1) then nil
else if degr(u3,var)=0
then rfactor := numr rfactor ./ multf(u3,denr rfactor)
else y := (multsq(rfactor,remsq(u2,!*f2q u3,var)) . car x) . y;
x := nil;
for each j in y do
if cddr j =1 then x := j . x
else x := append(pfpower(car j,cadr j,cddr j,var),x);
x := for each j in x
collect list('quotient,prepsq!* car j,
if cddr j=1 then prepf cadr j
else list('expt,prepf cadr j,cddr j));
if polypart then x := prepsq!* polypart . x;
setkorder xx;
return 'list . x
end;
% 11 May 98.
symbolic procedure simplimit u;
begin scalar fn,exprn,var,val,old,v,!*protfg;
if length u neq 4
then rerror(limit,1,
"Improper number of arguments to limit operator");
fn:= car u; exprn := cadr u; var := !*a2k caddr u; val := cadddr u;
!*protfg := t;
old := get('cot,'opmtch);
put('cot,'opmtch,
'(((!~x) (nil . t) (quotient (cos !~x) (sin !~x)) nil)));
v := errorset!*({'apply,mkquote fn,mkquote {exprn,var,val}},nil);
put('cot,'opmtch,old);
!*protfg := nil;
return if errorp v or (v := car v) = aeval 'failed then mksq(u,1)
else simp!* v
end;
endpatch;
% Numeric declarations.
fluid '(!*noequiv);
global '(iterations!*);
symbolic macro procedure dmx!: u;
subla('((plus2 . !:plus)(times2 . !:times)
(quotient . !:!:quotient)
(difference . !:difference)
(minus . !:minus)
(lessp . (lambda(a b)(!:minusp (!:difference a b))))
(greaterp . (lambda(a b)(!:minusp (!:difference b a))))
(abs . absf)),
cadr u);
patch numeric;
% 9 Aug 97.
symbolic procedure intrduni(e,p,acc);
dmx!:
begin scalar lo,hi,x,u,eps,i,il,int;
integer n,k;
x := car p; p:= cdr p;
lo := cadr x; hi := cddr x; x := car x;
lo:=force!-to!-dm lo; hi:=force!-to!-dm hi;
if hi=lo then return force!-to!-dm nil;
il := list intrdmkinterval(e,x,
(lo.intevaluate1(e,x,lo)),
(hi.intevaluate1(e,x,hi)),1);
int:=car lastpair car il;
loop:
k:=add1 k; n:=add1 n;
if remainder(n,25)=0 then intrdmess(n,il);
i:=car il; il:=cdr il; u:=intrdintersect(e,x,i);
if null u then
<<il:=i.il;
intrdabortmsg(list car cadr i,list x,intcollecteps il);
goto ready>>;
for each q in u do il := intrdsortin(q,il);
if k<5 then goto loop;
int:=intcollectint il;
k:=0; eps := intcollecteps il;
if eps>(acc* abs int) then goto loop;
ready:
return intcollectint il;
end;
symbolic procedure intrdmulti(e,p,acc);
dmx!:
begin scalar vars,u,eps,i,il,int;
integer n,k,dim;
if smemq('infinity,p) then
rederr "no support for infinity in multivariate integrals";
dim:=length p;
il:=intrdmkinitcube(e,p);
vars := car il; il:= cdr il;
loop:
k:=add1 k; n:=add1 n;
if remainder(n,25)=0 then intrdmess(n,il);
i:=car il; il:=cdr il; u:=intrdrefinecube(e,vars,i);
if null u then
<<il:=i.il;
intrdabortmsg(car cadr i,vars,intcollecteps il);
goto ready>>;
for each q in u do il := intrdsortin(q,il);
if k<5 then goto loop;
int:=intcollectint il;
k:=0; eps := intcollecteps il;
if eps>(acc* abs int) then goto loop;
ready:
return intcollectint il;
end;
symbolic procedure rdmineval u;
begin scalar e,vars,x,y,z,oldmode,p,!*noequiv;
oldmode:=steepdecedmode(nil,nil);
u := for each x in u collect reval x;
u := accuracycontrol(u,6,40);
e :=car u; u :=cdr u;
if eqcar(car u,'list) then
u := for each x in cdar u collect reval x;
for each x in u do
<<if eqcar(x,'equal) then
<<z:=cadr x; y:=caddr x>> else <<z:=x; y:=random 100>>;
vars:=z . vars; p := y . p>>;
vars := reversip vars; p := reversip p;
x := steepdeceval1(e,vars,p,'num_min);
steepdecedmode(t,oldmode);
return list('list, car x, 'list.
for each p in pair(vars,cdr x) collect
list('equal,car p,cdr p));
end;
symbolic procedure steepdec2(f,grad,vars,acc,x,mode);
dmx!:
begin scalar e0,e00,e1,e2,a,a1,a1a1,a2,a2a2,x1,x2,g,h,dx,
delta,limit,gold,gnorm,goldnorm,multi;
integer count,k,n;
n:=length grad; multi:=n>1; n:=10*n;
e00 := e0 := evaluate(f,vars,x);
a1 := 1;
init:
k:=0;
g := for each v in grad collect -evaluate(v,vars,x);
gnorm := normlist g; h:=g;
loop:
count := add1 count; k:=add1 k;
if count>iterations!* then
<<lprim "requested accuracy not reached within iteration limit";
goto ready>>;
a2 := nil;
l1:
x1 := list!+list(x, scal!*list(a1,h));
e1 := evaluate(f,vars,x1);
if e1>e0 then <<a2:=a1; x2:=x1; e2:=e1;
a1 := a1/2; goto l1>>;
if a2 then goto alph;
l2:
a2 := a1+a1;
x2 := list!+list(x, scal!*list(a2,h));
e2 := evaluate(f,vars,x2);
if e2<e1 then <<a1:=a2; e1:=e2; goto l2>>;
alph:
if abs(e1-e2)<acc then goto ready;
a1a1:=a1*a1; a2a2:=a2*a2;
a:= (a1a1-a2a2)*e0 + a2a2*e1 - a1a1*e2 ;
a:= a/((a1-a2)*e0 + a2*e1-a1*e2);
a:= a/2;
dx:=scal!*list(a,h); x:=list!+list(x, dx);
e0 := evaluate(f,vars,x);
if e0>e1 then
<<dx:=scal!*list(a1-a ,h); x:=list!+list(x,dx);
e0 := e1; dx:=scal!*list(a1,h)>>;
delta:= normlist dx;
steepdecprintpoint(count,x,delta,e0,gnorm);
update!-precision list(delta,e0,gnorm);
if k>n then goto init;
gold := g; goldnorm:= gnorm;
g := for each v in grad collect -evaluate(v,vars,x);
gnorm := normlist g;
if gnorm<limit then goto ready;
h := if not multi then g else
list!+list(g,scal!*list(gnorm/goldnorm,h));
if mode='num_min and gnorm > acc or
mode='root and e0 > acc then goto loop;
ready:
if mode='root and not(abs e0 < acc) then
<<lprim "probably fallen into local minimum of f^2";
return nil>>;
return e0 . x;
end;
endpatch;
% Ncpoly declarations.
fluid '(ncpi!-brackets!* ncpi!-comm!-rules!* ncpi!-name!-rules!*
ncpi!-names!* !*ncg!-right);
flag('(ncpi!-brackets!* ncpi!-comm!-rules!* ncpi!-name!-rules!*),
'share);
patch ncpoly;
% 15 Jun 96.
symbolic procedure ncpi!-setup u;
begin scalar vars,al,b,b0,f,m,rs,rn,na,rh,lh,s,x,y,w,!*evallhseqp;
if (w:=member('left,u)) or (w:=member('right,u)) then
<<u:=delete(car w,u); !*ncg!-right:=car w='right>>;
if length u > 2 then rederr "illegal number of arguments";
if ncpi!-name!-rules!* then
algebraic clearrules ncpi!-comm!-rules!*,ncpi!-name!-rules!*;
u:=sublis(ncpi!-names!*,u);
for each x in cdr listeval(car u,nil) collect
<<x:=reval x; y:={'nc!*,mkid('!_,x)}; na:=(y.x).na;
rn:={'replaceby,x,y}.rn; al:=(x.y).al;vars:=append(vars,{y})>>;
ncpi!-names!*:=na;
ncpi!-name!-rules!*:='list.rn;
m:=for i:=1:length vars -1 join
for j:=i+1:length vars collect nth(vars,i).nth(vars,j);
if cdr u then ncpi!-brackets!*:=listeval(cadr u,nil);
if null ncpi!-brackets!* then
rederr "commutator relations missing";
for each b in cdr ncpi!-brackets!* do
<<b0:=sublis(al,b);
w:= eqcar(b0,'equal) and (lh:=cadr b0) and (rh:=caddr b0)
and eqcar(lh,'difference) and (f:=cadr lh) and (s:=caddr lh)
and eqcar(f,'times) and eqcar(s,'times)
and (x:=cadr f) and (y:=caddr f) and member(x,vars)
and member(y,vars)
and x=caddr s and y=cadr s;
if not w then typerr(b,"commutator relation");
if member(x.y,m) then <<w:=x;x:=y;y:=w; rh:={'minus,rh}>>;
m:=delete(y.x,m);
rs:={'replaceby,{'times,x,y},{'plus,{'times,y,x},rh}} .rs>>;
ncdsetup!*{'list.vars,'list.rs};
apply('korder,{vars});
apply('order,{vars});
for each c in m do
rs:={'replaceby,{'times,cdr c,car c},{'times,car c,cdr c}}. rs;
ncpi!-comm!-rules!*:='list.rs;
algebraic let ncpi!-comm!-rules!*,ncpi!-name!-rules!*;
end;
endpatch;
% Poly declarations.
fluid '(!*algint !*!*processed);
patch poly;
% 27 Dec 95, 28 Mar 96.
symbolic procedure repartf u;
(if domainp u
then if atom u then u
else if get(car u,'cmpxfn)
then int!-equiv!-chk(car u . cadr u .
cadr apply1(get(car u,'i2d),0))
else u
else if mvar u eq 'i then repartf red u
else addf(multpf(lpow u,repartf lc u),repartf red u))
where u = reorder u where kord!* = 'i . kord!*;
% 27 Dec 95.
symbolic procedure impartf u;
(if domainp u
then if atom u then nil
else if get(car u,'cmpxfn)
then int!-equiv!-chk(car u . cddr u .
cadr apply1(get(car u,'i2d),0))
else nil
else if mvar u eq 'i then addf(lc u,impartf red u)
else if null dmode!* then impartf red u
else addf(multpf(lpow u,impartf lc u),impartf red u))
where u = reorder u where kord!* = 'i . kord!*;
% 6 Jan 96.
symbolic procedure factorf u;
(begin scalar m!-image!-variable,new!-korder,old!-korder,sign,v,w;
if domainp u then return list u;
new!-korder:=kernord(u,nil);
if !*kernreverse then new!-korder:=reverse new!-korder;
old!-korder:=setkorder new!-korder;
u := reorder u;
if minusf u then <<sign := not sign; u := negf u>>;
v := comfac u;
u := quotf1(u,cdr v);
if domainp u then errach list("Improper factors in factorf");
w := car v;
v := fctrf1 cdr v;
if w then v := car v . (!*k2f car w . cdr w) . cdr v;
m!-image!-variable := mvar u;
u :=
distribute!.multiplicity(factorize!-primitive!-polynomial u,1);
setkorder old!-korder;
if sign then u := (negf caar u . cdar u) . cdr u;
u := fac!-merge(v,1 . u);
return car u . for each w in cdr u collect (reorder car w . cdr w)
end) where current!-modulus = current!-modulus;
% 11 Jan 96.
symbolic procedure diffp(u,v);
begin scalar n,w,x,y,z; integer m;
n := cdr u;
u := car u;
if u eq v and (w := 1 ./ 1) then go to e
else if atom u then go to f
else if (not atom car u and (w:= difff(u,v)))
or (car u eq '!*sq and (w:= diffsq(cadr u,v)))
then go to c
else if x := get(car u,'dfform) then return apply3(x,u,v,n)
else if x:= get(car u,dfn_prop u) then nil
else if car u eq 'plus and (w := diffsq(simp u,v))
then go to c
else go to h;
y := x;
z := cdr u;
a: w := diffsq(simp car z,v) . w;
if caar w and null car y then go to h;
y := cdr y;
z := cdr z;
if z and y then go to a
else if z or y then go to h;
y := reverse w;
z := cdr u;
w := nil ./ 1;
b: if caar y
then w := addsq(multsq(car y,simp subla(pair(caar x,z),
cdar x)),
w);
x := cdr x;
y := cdr y;
if y then go to b;
c:
e: if (x := atsoc(u,wtl!*))
then w := multpq('k!* .** (-cdr x),w);
m := n-1;
return rationalizesq if n=1 then w
else if flagp(dmode!*,'convert)
and null(n := int!-equiv!-chk
apply1(get(dmode!*,'i2d),n))
then nil ./ 1
else multsq(!*t2q((u .** m) .* n),w);
f: if not depends(u,v)
and (not (x:= atsoc(u,powlis!*))
or not depends(cadddr x,v))
then return nil ./ 1;
w := list('df,u,v);
w := if x := opmtch w then simp x else mksq(w,1);
go to e;
h: if car u eq 'df
then if depends(cadr u,v)
and not(cadr u eq v and not depends(v,caddr u))
then
<<if !*fjwflag and eqcar(cadr u, 'int) then
if caddr cadr u eq v then
<< w := 'df . cadr cadr u . cddr u; go to j >>
else if !*allowdfint and
not_df_p(w := diffsq(simp!* cadr cadr u, v))
then <<
w := aeval{'int, mk!*sq w, caddr cadr u} . cddr u;
go to j >>;
if (x := find_sub_df(w:= cadr u . derad(v,cddr u),
get('df,'kvalue)))
then <<w := simp car x;
for each el in cdr x do
for i := 1:cdr el do
w := diffsq(w,car el);
go to e>>
else w := 'df . w
>>
else return nil ./ 1
else w := list('df,u,v);
j: if (x := opmtch w) then w := simp x
else if not depends(u,v) then return nil ./ 1
else w := mksq(w,1);
go to e
end;
% 17 Jan 96.
symbolic procedure exptsq(u,n);
begin scalar x;
if n=1 then return u
else if n=0
then return if null numr u then rerror(poly,4," 0**0 formed")
else 1 ./ 1
else if null numr u then return u
else if n<0 then return simpexpt list(mk!*sq u,n)
else if null !*exp
then return mksfpf(numr u,n) ./ mksfpf(denr u,n)
else if kernp u then return mksq(mvar numr u,n)
else if denr u=1 then return exptf(numr u,n) ./ 1
else if domainp numr u
then x := multsq(!:expt(numr u,n) ./ 1,1 ./ exptf(denr u,n))
else <<x := u;
while (n := n-1)>0
do x := multf(numr u,numr x) ./ multf(denr u,denr x);
x := canonsq x>>;
if null cdr x then rerror(poly,101,"Zero divisor");
return x
end;
% 11 Feb 96.
symbolic procedure dcombine(u,v,fn);
<<u := if atom u
then apply2(get(car v,fn),apply1(get(car v,'i2d),u),v)
else if atom v
then apply2(get(car u,fn),u,apply1(get(car u,'i2d),v))
else if car u eq car v then apply2(get(car u,fn),u,v) else
(<<if (not(x and atom x)
or (get(du,'cmpxfn) and not get(dv,'cmpxfn)
or du memq dl and not(dv memq dl)) and dv neq '!:ps!:)
and not flagp(dv,'noconvert) then
<<if du memq dml and dv eq '!:rd!:
or du eq '!:rd!: and dv memq dml then
<<u := apply1(get(du,'!:cr!:),u); du := '!:cr!:>>
else if du eq '!:rn!: and dv eq '!:gi!:
or du eq '!:gi!: and dv eq '!:rn!: then
<<u := apply1(get(du,'!:crn!:),u); du := '!:crn!:>>;
v := apply1(get(dv,du),v); x := get(du,fn)>>
else <<u := apply1(x,u); x := get(dv,fn)>>;
apply2(x,u,v)>>
where x=get(du,dv),dml='(!:crn!: !:gi!:),dl='(!:rd!: !:cr!:))
where du=car u,dv=car v;
if !*rounded and !*roundall and not atom u then
int!-equiv!-chk
if (v := car u) eq '!:rn!: and cddr u neq 1 then !*rn2rd u
else if v eq '!:crn!: and (cdadr u neq 1 or cdddr u neq 1)
then !*crn2cr u
else u
else if fn eq 'divide then
int!-equiv!-chk car u . int!-equiv!-chk cdr u
else int!-equiv!-chk u>>;
% 6 Mar 96, 16 Jul 96, 17 Jan 98.
symbolic procedure gcdf(u,v);
begin scalar !*exp, !*rounded;
!*exp := t;
u := if domainp u or domainp v or not !*ezgcd
or dmode!*
or free!-powerp u or free!-powerp v
then gcdf1(u,v)
else ezgcdf(u,v);
return if minusf u then negf u else u
end;
symbolic procedure gcdnx(u,v);
if domainp u then gcdfd(u,v)
else if domainp v then gcdfd(v,u)
else if null red u then gcdnx(lc u,v)
else gcdnx(lc u,gcdnx(red u,v));
% 21 May 96.
symbolic procedure mksfpf(u,n);
(if n=1 then x
else if domainp x then !:expt(x,n)
else if n>=0 and onep lc x and null red x
then (((if z and cdr z<=m then nil else !*p2f mksp(y,m))
where z=assoc(y,asymplis!*)) where m=ldeg x*n,y=mvar x)
else exptf2(x,n))
where x=mkprod u;
% 22 May 96, 20 Jul 98.
symbolic procedure multf(u,v);
begin scalar x,y;
a: if null u or null v then return nil
else if u=1 then return v
else if v=1 then return u
else if domainp u then return multd(u,v)
else if domainp v then return multd(v,u)
else if not(!*exp or ncmp!* or wtl!* or x)
then <<u := mkprod u; v := mkprod v; x := t; go to a>>;
x := mvar u;
y := mvar v;
if noncomfp v and (noncomp x or null !*!*processed)
then return multfnc(u,v)
else if x eq y
then <<if not fixp ldeg u or not fixp ldeg v
then x := x .** reval list('plus,ldeg u,ldeg v)
else x := mkspm(x,ldeg u+ldeg v);
y := addf(multf(red u,v),multf(!*t2f lt u,red v));
return if null x or null(u := multf(lc u,lc v))
then <<!*asymp!* := t; y>>
else if x=1 then addf(u,y)
else if null !*mcd then addf(!*t2f(x .* u),y)
else x .* u .+ y>>
else if ordop(x,y)
then <<x := multf(lc u,v);
y := multf(red u,v);
return if null x then y else lpow u .* x .+ y>>;
x := multf(u,lc v);
y := multf(u,red v);
return if null x then y else lpow v .* x .+ y
end;
symbolic procedure noncomfp u;
not domainp u and (noncomp mvar u or noncomfp lc u or noncomfp red u);
symbolic procedure multfnc(u,v);
begin scalar x,y;
x := multf(lc u,!*t2f lt v);
if null x then nil
else if not domainp x and mvar x eq mvar u
then x := addf(if null (y := mkspm(mvar u,ldeg u+ldeg x))
then nil
else if y = 1 then lc x
else !*t2f(y .* lc x),
multf(!*p2f lpow u,red x))
else if noncomp mvar u then x := !*t2f(lpow u .* x)
else x := multf(!*p2f lpow u,x) where !*!*processed=t;
return addf(x,addf(multf(red u,v),multf(!*t2f lt u,red v)))
end;
% 26 May 96, 1 Sep 97, 31 Jan 98.
symbolic procedure rationalizesq u;
begin scalar !*structure,!*sub2,v,x;
if x := get(dmode!*,'rationalizefn) then u := apply1(x,u);
v := subs2q u;
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;
v := subs2q v;
if not domainp numr quotsq(u,v)
then rationalizesq v else v>>
else u
end;
% 21 Jul 96.
symbolic procedure mkprod u;
begin scalar w,x,y,z,!*exp,!*sub2;
if null u or kernlp u then return u;
!*sub2 := t;
x := subs2(u ./ 1);
if denr x neq 1 then return u
else if numr x neq u
then <<u := numr x; if null u or kernlp u then return u>>;
!*exp := t;
w := ckrn u;
u := quotf(u,w);
x := expnd u;
if null x or kernlp x then return multf(w,x);
if !*mcd and (!*sqfree or !*factor or !*gcd)
then y := fctrf x
else <<y := ckrn x; x := quotf(x,y); y := list(y,x . 1)>>;
if cdadr y>1 or cddr y
then <<z := car y;
for each j in cdr y do
z := multf(mksp!*(car j,cdr j),z)>>
else if not !*group and tmsf u>tmsf caadr y
then z := multf(mksp!*(caadr y,cdadr y),car y)
else z := mksp!*(u,1);
return multf(w,z)
end;
% 27 Apr 96, 27 Jul 96.
symbolic procedure gcdf2(u,v);
begin scalar asymplis!*,w,z;
if not num!-exponents u or not num!-exponents v then return 1;
if !*gcd and length(w := kernord(u,v))>1
then <<w := list setkorder w;
u := reorder u;
v := reorder v>>
else w := nil;
if mvar u eq mvar v
then begin scalar x,y;
x := comfac u;
y := comfac v;
z := gcdf1(cdr x,cdr y);
u := quotf1(u,comfac!-to!-poly x);
v := quotf1(v,comfac!-to!-poly y);
if !*gcd then z := multf(gcdk(u,v),z)
else if v and quotf1(u,v) then z := multf(v,z)
else if u and quotf1(v,u) then z := multf(u,z);
if car x and car y
then if pdeg car x>pdeg car y
then z := multpf(car y,z)
else z := multpf(car x,z)
end
else if ordop(mvar u,mvar v) then z := gcdf1(cdr comfac u,v)
else z := gcdf1(cdr comfac v,u);
if w then <<setkorder car w; z := reorder z>>;
return z
end;
% 27 Jul 96.
symbolic procedure quotf1(p,q);
if null p then nil
else if p=q then 1
else if q=1 then p
else if domainp q then quotfd(p,q)
else if domainp p then nil
else if mvar p eq mvar q
then begin scalar u,v,w,x,xx,y,z,z1; integer n;
a:if idp(u := rank p) or idp(v := rank q) or u<v then return nil;
u := lt!* p;
v := lt!* q;
w := mvar q;
x := quotf1(tc u,tc v);
if null x then return nil;
n := tdeg u-tdeg v;
if n neq 0 then y := w .** n;
xx := multf(negf x,red q);
p := addf(red p,if n=0 then xx else multpf(y,xx));
if p and (domainp p or not(mvar p eq w)) then return nil
else if n=0 then go to b;
z := aconc!*(z,y .* x);
if null p then return if z1 then rnconc(z,z1) else z;
go to a;
b: if null p then return rnconc(z,x)
else if !*mcd then return nil
else z1 := x;
go to a
end
else if ordop(mvar p,mvar q) then quotk(p,q)
else nil;
% 4 Aug 96;
symbolic procedure rnconc(u,v);
if null u then v
else if noncomfp u and noncomfp v then multf(u,v)
else begin scalar w;
w := u;
while cdr w do <<w := cdr w>>;
rplacd(w,v);
return u
end;
symbolic procedure noncomfp u;
not domainp u and (noncomp mvar u or noncomfp lc u or noncomfp red u);
% 11 Sep 96.
symbolic procedure setcmpxmode(u,bool);
begin scalar x,y;
x := get(u,'tag);
if u eq 'complex
then if null dmode!*
then return if null bool then nil
else <<put('i,'idvalfn,'mkdgi);
setdmode1('complex,bool)>>
else if null bool
then return if null !*complex then nil
else if get(dmode!*,'dname) eq 'complex
then <<remprop('i,'idvalfn);
setdmode1('complex,nil)>>
else <<remprop('i,'idvalfn);
setdmode1(get(get(dmode!*,'realtype),'dname),
t)>>
else if dmode!* eq '!:gi!: or get(dmode!*,'realtype)
then return nil
else if not (y := get(dmode!*,'cmpxtype))
then dmoderr(dmode!*,x)
else <<put('i,'idvalfn,get(y,'ivalue));
return setdmode1(get(y,'dname),bool)>>
else if null bool then
if u eq (y := get(get(dmode!*,'realtype),'dname)) then
<<put('i,'idvalfn,'mkdgi); return setdmode1('complex,t)>>
else if null y then return nil
else offmoderr(u,y)
else <<u := get(u,'tag);
y := get(u,'cmpxtype);
if null y then dmoderr(u,'!:gi!:);
put('i,'idvalfn,get(y,'ivalue));
return setdmode1(get(y,'dname),bool)>>
end;
% 3 Nov 96.
symbolic procedure lcofeval u;
begin scalar kern,x,y;
if null u or null cdr u or not null cddr u
then rerror(poly,280,
"LCOF called with wrong number of arguments");
kern := !*a2k cadr u;
u := simp!* car u;
y := denr u;
tstpolyarg(y,u);
u := numr u;
if domainp u then return if null u then 0 else mk!*sq (u . 1)
else if mvar u eq kern then return !*ff2a(lc u,y);
x := updkorder kern;
u := reorder u;
if mvar u eq kern then u := lc u;
setkorder x;
return if null u then 0 else !*ff2a(u,y)
end;
% 22 Dec 96.
symbolic procedure simpdf u;
begin scalar v,x,y,z;
if null subfg!* then return mksq('df . u,1);
v := cdr u;
u := simp!* car u;
a: if null v or null numr u then return u;
x := if null y or y=0 then simp!* car v else y;
if denr x neq 1 or atom numr x
then typerr(prepsq x,"kernel or integer")
else (if domainp z
then if get(car z,'domain!-diff!-fn)
then begin scalar dmode!*,alglist!*;
x := prepf z;
if null prekernp x
then typerr(x,"kernel")
end
else typerr(prepf z,"kernel")
else if null red z and lc z = 1 and ldeg z = 1
then x := mvar z
else typerr(prepf z,"kernel")) where z = numr x;
v := cdr v;
if null v then <<u := diffsq(u,x); go to a>>;
y := simp!* car v;
if null numr y then <<v := cdr v; y := nil; go to a>>
else if not(z := d2int y) then <<u := diffsq(u,x); go to a>>;
v := cdr v;
for i:=1:z do u := diffsq(u,x);
y := nil;
go to a
end;
symbolic procedure d2int u;
if denr u neq 1 then nil
else if numberp(u := numr u) then u
else if not domainp u or not(car u eq '!:rd!:) then nil
else (if abs(float x - u)<!!fleps1 then x else nil)
where x=fix u where u=rd2fl u;
% 30 Jan 97.
symbolic procedure !:quotient(u,v);
if null u or u=0 then nil
else if null v or v=0 then rerror(poly,12,"Zero divisor")
else if atom u and atom v
% We might also check that remainder is zero in integer case.
then if null dmode!* then u/v else
(if atom recipv then u*recipv else dcombine(u,recipv,'times))
where recipv=!:recip v
else dcombine(u,v,'quotient);
% 9 Apr 97.
symbolic procedure comfac p;
begin scalar x,y;
if flagp(dmode!*,'field) and ((x := lnc p) neq 1)
then p := multd(!:recip x,p);
if null red p then return lt p;
x := lc p;
y := mvar p;
a: p := red p;
if degr(p,y)=0
then return nil
. if domainp p or not(noncomp y and noncomp mvar p)
then gcdf(x,p)
else 1
else if null red p then return lpow p . gcdf(x,lc p)
else x := gcdf(lc p,x);
go to a
end;
% 24 Dec 97, 18 Feb 98, 31 Mar 98.
symbolic procedure sqfrf u;
begin integer n; scalar !*gcd,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)>>
else if (w := get(dmode!*,'units)) and (w := assoc(y:= lnc u,w))
then <<u := multd(cdr w,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 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;
% 18 Feb 98.
symbolic procedure factor!-prim!-f u;
begin scalar v,w,x,y;
if ncmp!* then return list(1,u . 1);
if dmode!* and (x := get(dmode!*,'sqfrfactorfn))
then if !*factor then v := apply1(x,u) else v := list(1,u . 1)
else if flagp(dmode!*,'field) and ((w := lnc u) neq 1)
then v := w . sqfrf multd(!:recip w,u)
else if (w := get(dmode!*,'units)) and (w := assoc(y := lnc u,w))
then v := y . sqfrf multd(cdr w,u)
else v := 1 . sqfrf u;
if x and (x eq get(dmode!*,'factorfn))
then return v;
w := list car v;
for each x in cdr v
do w := fac!-merge(factor!-prim!-sqfree!-f x,w);
return w
end;
endpatch;
patch rlisp;
% 4 Dec 95.
symbolic procedure yesp1;
begin scalar bool,x,y;
a: x := readch();
if x eq !$eol!$ then go to a
else if x eq !$eof!$ then eval '(bye)
else if (y := x memq '(!y !Y)) or x memq '(!n !N)
then return y
else if null bool then <<prin2t "Type Y or N"; bool := t>>;
go to a
end;
% 18 Dec 96.
if 'csl memq lispsystem!*
then <<put('gc,'simpfg,'((t (verbos t)) (nil (verbos nil))));
switch gc>>;
% 18 Feb 96.
symbolic procedure getsetvarlis u;
if null u then nil
else if atom u then errach list("getsetvarlis",u)
else if atom car u then car u . getsetvarlis cdr u
else if caar u memq '(setel setk)
then getsetvarlis cadar u . getsetvarlis cdr u
else if caar u eq 'setq then mkquote cadar u . getsetvarlis cdr u
else car u . getsetvarlis cdr u;
symbolic procedure getsetvars u;
if atom u then nil
else if car u memq '(setel setk)
then getsetvarlis cadr u . getsetvars caddr u
else if car u eq 'setq then mkquote cadr u . getsetvars caddr u
else nil;
% 31 Dec 97.
symbolic procedure !*!*a2s(u,vars);
if null u or constantp u and null fixp u
or intexprnp(u,vars) and null !*composites
and null current!-modulus
or flagpcar(u,'nochange) and not(car u eq 'getel)
then u
else if smember('random,u) then
list(list('lambda,'(!*uncached),list(!*!*a2sfn,u)),t)
else list(!*!*a2sfn,u);
% 1 Jan 98.
symbolic procedure formc!*(u,vars,mode);
begin scalar !*!*a2sfn;
!*!*a2sfn := 'revalx;
return formc(u,vars,mode)
end;
symbolic procedure revalx u;
reval if not atom u and not atom car u then prepf u else u;
endpatch;
patch scope;
% 17 Feb 96.
symbolic procedure primat;
begin scalar l;
l := linelength(120);
terpri();
prin2 "Sumscheme :";
prischeme('plus);
terpri();
terpri();
terpri();
prin2 "Productscheme :";
prischeme('times);
linelength(l);
end;
endpatch;
% Solve declarations.
fluid '(!*fullroots !*solvealgp !*solvesingular !!gcd vars!*
!*!*norootvarrenamep!*!*);
global '(assumptions);
patch solve;
% 20 Apr 96.
symbolic procedure solvemixedsys(exlis,varlis);
if null cadr(exlis := siftnonlnr(exlis,varlis)) then
solvelnrsys(car exlis,varlis)
else if null car exlis then
solvenonlnrsys(cadr exlis,varlis)
else
begin scalar x,y,z;
x := solvelnrsys(car exlis,varlis) where !*arbvars = nil;
if car x = 'inconsistent then return x
else x := cadr x;
z := pair(cadr x,foreach ex in car x collect mk!*sq ex);
exlis := foreach ex in cadr exlis join
if ex := numr subs2 resimp subf(ex,z) then {ex};
varlis := setdiff(varlis,cadr x);
y := solvemixedsys(exlis,varlis);
if car y memq {'inconsistent,'singular,'failed,nil} then return y
else return t . foreach s in cdr y collect
<<z := foreach pr in pair(cadr s,car s) join
if not smemq('root_of,cdr pr) then
{car pr . mk!*sq cdr pr};
{append(car s,foreach ex in car x collect subsq(ex,z)),
append(cadr s,cadr x),caddr s}>>;
end;
% 12 Dec 95, 9 May 96.
symbolic procedure solvesq1 (ex,var,mul);
begin scalar e1,oldkorder,x1,y,z; integer mu;
ex := numr ex;
if null(x1 := topkern(ex,var)) then return nil;
oldkorder := setkorder list var;
e1 := reorder ex;
setkorder oldkorder;
if !*modular then
<<load!_package 'modsr; return msolvesys({e1},x1,nil)>>;
if mvar e1 = var and null cdr x1 and ldeg e1 =1
then return {{{quotsq(negf reorder red e1 ./1,
reorder lc e1 ./ 1)},
{var},mul}};
ex := if !*rounded then {1,ex . 1} else fctrf ex;
if domainp car ex then ex := cdr ex
else ex := (car ex . 1) . cdr ex;
for each j in ex do
<<e1 := car j;
x1 := topkern(e1,var);
mu := mul*cdr j;
if null !*rounded and length x1=1 and length kernels e1<5
and length(y := decomposef1(e1,nil))>1
and (y := solvedecomp(reverse y,car x1,mu))
then z := append(y,z)
else if (degr(y := reorder e1,var) where kord!*={var}) = 1
and not smember(var,delete(var,x1))
then <<y := {{quotsq(!*f2q negf reorder red y,
!*f2q reorder lc y)},
{var},mu};
z := y . z>>
else if x1
then z := solnsmerge(
if null cdr x1 then solve1(e1,car x1,var,mu)
else if (y := principle!-of!-powers!-soln(e1,x1,var,mu))
neq 'unsolved
then y
else if not smemq('sol,x1 := solve!-apply!-rules(e1,var))
then solvesq(x1,var,mu)
else mkrootsof(e1 ./ 1,var,mu),
z)>>;
return z
end;
symbolic procedure solve11(e1,x1,var,mu);
begin scalar !*numval,b,coefs,hipow,w; integer n;
if null !*fullroots and null !*rounded and numrdeg(e1,var)>2
then return mkrootsof(e1 ./ 1,var,mu);
!*numval := t;
coefs:= errorset!*(list('solvecoeff,mkquote e1,mkquote x1),nil);
if atom coefs or atom x1 and x1 neq var
then return mkrootsof(e1 ./ 1,var,mu);
coefs := car coefs;
n:= !!gcd;
hipow := car(w:=car reverse coefs);
if not domainp numr cdr w then
assumptions:=append(assumptions,{prepf numr cdr w});
if not domainp denr cdr w then
assumptions:=append(assumptions,{prepf denr cdr w});
if hipow = 1
then return begin scalar lincoeff,y,z;
if null cdr coefs then b := 0
else b := prepsq quotsq(negsq cdar coefs,cdadr coefs);
if n neq 1 then b := list('expt,b,list('quotient,1,n));
for k := 0:n-1 do
<<lincoeff := list('times,b,
mkexp list('quotient,list('times,k,2,'pi),n));
lincoeff := simp!* lincoeff where alglist!* = nil . nil;
if x1=var
then y := solnmerge(list lincoeff,list var,mu,y)
else if not idp(z := car x1)
then typerr(z,"solve operator")
else if z := get(z,'solvefn)
then y := append(apply1(z,list(cdr x1,var,mu,lincoeff))
,y)
else if (z := get(car x1,'inverse))
then y := append(solvesq(subtrsq(simp!* cadr x1,
simp!* list(z,mk!*sq lincoeff)),
var,mu),y)
else y := list(list subtrsq(simp!* x1,lincoeff),nil,mu)
. y>>;
return y
end
else if hipow=2
then return <<x1 := exptsq(simp!* x1,n);
w := nil;
for each j in solvequadratic(getcoeff(coefs,2),
getcoeff(coefs,1),getcoeff(coefs,0))
do w :=
solnsmerge(solvesq(subtrsq(x1,j),var,mu),w);
w>>
else return solvehipow(e1,x1,var,mu,coefs,hipow)
end;
symbolic procedure solnsmerge(u,v);
if null u then v
else solnsmerge(cdr u,solnmerge(caar u,cadar u,caddar u,v));
symbolic procedure solvehipow(e1,x1,var,mu,coefs,hipow);
begin scalar b,c,d,f,rcoeffs,z;
f:=(hipow+1)/2;
d:=exptsq(simp!* x1,!!gcd);
rcoeffs := reverse coefs;
return if solve1test1(coefs,rcoeffs,f)
then if f+f=hipow+1
then <<c:=addsq(d, 1 ./ 1);
append(solvesq(c,var,mu),
solvesq(quotsq(e1 ./ 1, c),var,mu))>>
else <<coefs := putcoeff(coefs,0,2 ./ 1);
coefs := putcoeff(coefs,1,simp!* '!!x);
c:=addsq(multsq(getcoeff(coefs,f+1),
getcoeff(coefs,1)),
getcoeff(coefs,f));
for j:=2:f do <<
coefs := putcoeff(coefs,j,
subtrsq(multsq(getcoeff(coefs,1),
getcoeff(coefs,j-1)),
getcoeff(coefs,j-2)));
c:=addsq(c,multsq(getcoeff(coefs,j),
getcoeff(coefs,f+j)))>>;
for each j in solvesq(c,'!!x,mu) do z :=
solnsmerge(
solvesq(addsq(1 ./ 1,multsq(d,subtrsq(d,caar j))),
var,caddr j),z);
z>>
else if solve1test2(coefs,rcoeffs,f)
then <<c:=addsq(d,(-1 ./1));
b := solvesq(c,var,mu);
e1 := quotsq(e1 ./ 1, c);
if f+f = hipow
then <<c := addsq(d,(1 ./ 1));
b := append(solvesq(c,var,mu),b);
e1 := quotsq(e1,c)>>;
append(solvesq(e1,var,mu),b)>>
else if !*rounded and univariatep e1
then reversip solveroots(e1,var,mu)
else if null !*fullroots then mkrootsof(e1 ./ 1,var,mu)
else if hipow=3
then <<for each j in solvecubic(getcoeff(coefs,3),
getcoeff(coefs,2),
getcoeff(coefs,1),
getcoeff(coefs,0))
do z := solnsmerge(solvesq(subtrsq(d,j),var,mu),z);
z>>
else if hipow=4
then <<for each j in solvequartic(getcoeff(coefs,4),
getcoeff(coefs,3),
getcoeff(coefs,2),
getcoeff(coefs,1),
getcoeff(coefs,0))
do z := solnsmerge(solvesq(subtrsq(d,j),var,mu),z);
z>>
else mkrootsof(e1 ./ 1,var,mu)
end;
% 28 Dec 95.
symbolic procedure all_real(a,b,c);
begin scalar dmode,!*ezgcd,!*msg,!*numval;
dmode := dmode!*;
!*numval := t;
on complex,rounded;
a := real_1(a := resimp a) and real_1(b := resimp b)
and real_1(c := resimp c);
off rounded,complex;
if dmode then onoff(get(dmode,'dname),t);
return a
end;
% 6 Jan 96.
symbolic procedure solveroots(ex,var,mu);
begin scalar x,y;
x := reval list('roots_at_prec,mk!*sq(ex ./ 1));
if not(car x eq 'list)
then errach list("incorrect root format",ex);
for each z in cdr x do
y := solnsmerge(
solvesq(if not(car z eq 'equal)
then errach list("incorrect root format",ex)
else simpplus {cadr z,{'minus,caddr z}},
var,mu),
y);
return y
end;
% 12 Jan 96, 4 Jun 96.
symbolic procedure solvenonlnrsys(sys,uv);
begin scalar q,r,s,tag,!*expandexpt;
s := '(nil);
if solve!-psysp(sys,uv) then s := {sys} else
for each p in sys do
<<r := nil;
for each q in cdr fctrf p do
if topkernlis(car q,uv)
then for each u in s do r := (car q . u) . r;
s := r>>;
tag := 'failed; r := nil;
for each u in s do
<<u := solvenonlnrcollectexpt u;
q := solvenonlnrsys1(u,uv);
if eqcar(q,'failed) then q := solvenonlnrsyssep(u,uv);
if eqcar(q,'failed) then q := solvenonlnrsyslin(u,uv,nil);
if eqcar(q,'not) then q := solvenonlnrsyslin(u,uv,t);
if eqcar(q,'not) then q := '(failed);
if car q and car q neq 'failed then tag := car q;
q := if car q neq 'failed then cdr q else
for each j in u collect {{j ./ 1},nil,1};
r := union(q,r)>>;
return if tag eq 'inconsistent or tag eq 'failed then {tag}
else tag . r
end;
symbolic procedure topkernlis(u,v);
v and (topkern(u,car v) or topkernlis(u,cdr v));
% 26 Feb 96.
symbolic procedure solve0(elst,xlst);
begin scalar !*exp,!*notseparate,w; integer neqn;
!*exp := t;
elst := for each j in solveargchk elst collect simp!* !*eqn2a j;
neqn := length elst;
if null xlst
then <<vars!* := solvevars elst;
terpri();
if null vars!* then nil
else if cdr vars!*
then <<prin2!* "Unknowns: "; maprin('list . vars!*)>>
else <<prin2!* "Unknown: "; maprin car vars!*>>;
terpri!* nil>>
else <<xlst := solveargchk xlst;
vars!* := for each j in xlst collect !*a2k j>>;
if length vars!* = 0
then rerror(solve,3,"SOLVE called with no variables");
if neqn = 1 and length vars!* = 1
then if null numr car elst
then return if !*solvesingular
then {{{!*f2q makearbcomplex()},vars!*,1}}
else nil
else if solutionp(w := solvesq(car elst,car vars!*,1))
or null !*solvealgp
or univariatep numr car elst
then return w;
elst := for each j in elst collect numr j;
w := solvesys(elst,vars!*);
if null w then return nil;
if car w memq {'t,'inconsistent,'singular} then return cdr w
else if car w eq 'failed or null car w
then return for each j in elst collect list(list(j ./ 1),nil,1)
else errach list("Improper solve solution tag",car w)
end;
% 23 Apr 96.
switch varopt;
% 9 May 96.
symbolic procedure varsift(a,var);
if atom a then
if not(null a or numberp a or a eq var) then list a else nil
else if get(car a,'dname) then nil
else if car a memq '(arbint arbcomplex)
or (get(car a,'simpfn) eq 'simpiden and not smember(var,a))
then list a
else if car a eq '!*sq then varsift(prepsq cadr a,var)
else for each c in cdr a join varsift(c,var);
% 20 Jan 97.
symbolic procedure siftnonlnr(exlis,varlis);
begin scalar lin,nonlin;
foreach ex in exlis do
if ex then
if nonlnr(ex,varlis) then nonlin := ex . nonlin
else lin := ex . lin;
return {reversip lin,reversip nonlin};
end;
% 31 Jan 97.
symbolic procedure solve!-clean!-info(fl,flg);
begin scalar r,w,p;
for each form in cdr fl do if not p then
if constant_exprp(form:=reval form) then
(if not flg then p:= r:={1})
else
if flg then
for each w in cdr fctrf numr simp form do
<<w := absf car w;
if not member(w,r) then r := w . r>>
else
<<w:= absf numr simp{'nprimitive,form};
if not domainp w then w := sqfrf!* w;
for each z in r do if w then
if null cdr qremf(z,w) then r:=delete(z,r)
else if null cdr qremf(w,z) then w := nil;
if w then r := w . r>>;
return 'list . for each q in r collect prepf q
end;
% 20 Feb 97.
symbolic procedure safe!-modrecip u;
begin scalar q,!*msg,!*protfg;
!*msg := nil; !*protfg := t;
if eqcar(u,'!:mod!:) then u := cdr u;
q := errorset({'general!-modular!-reciprocal, u},nil,nil);
erfg!* := nil;
return if errorp q then nil else car q
end;
% 23 Feb 97.
symbolic procedure solvelnrsys(exlis,varlis);
begin scalar w,method;
if w := solvesparsecheck(exlis,varlis) then exlis := w
else exlis := exlis . varlis;
if null !*cramer and null exptexpflistp car exlis
then method := 'solvebareiss
else method := 'solvecramer;
exlis := apply2(method,car exlis,cdr exlis)
where bareiss!-step!-size!* = if w then 4 else 2;
return solvesyspost(exlis,varlis);
end;
% 29 Dec 97.
symbolic procedure mkrootsof(e1,var,mu);
begin scalar x,name;
x := if idp var then var else 'q_;
name := !*!*norootvarrenamep!*!* or mkrootsoftag();
if not !*!*norootvarrenamep!*!*
then while smember(x,e1) do
x := intern compress append(explode x,explode '!_);
e1 := prepsq!* e1;
if x neq var then e1 := subst(x,var,e1);
return list list(list !*kk2q list('root_of,e1,x,name),list var,mu)
end;
endpatch;
patch specfn;
% 10 Jan 96.
algebraic procedure compute_int_functions(x,f);
begin scalar pre,!*uncached,scale,term,n,precis,result,interm;
pre := precision 0; precision pre;
precis := 10^(-2 * pre);
lisp (!*uncached := t);
if f = Si then
if x < 0 then result :=
- compute_int_functions(-x,f) else
<< n:=1; term := x; result := x;
while abs(term) > precis do
<< term := -1 * (term * x*x)/(2n * (2n+1));
result := result + term/(2n+1);
n := n + 1>>; >>
else if f = Ci then
if x < 0 then result :=
- compute_int_functions(-x,f) -i*pi else
<< n:=1; term := 1; result := euler!*constant + log(x);
while abs(term) > precis do
<< term := -1 * (term * x*x)/((2n-1) * 2n);
result := result + term/(2n);
n := n + 1>>; >>
else if f = Ei then
<< n:=1; term := 1; result := euler!*constant + log(x);
while abs(term) > precis do
<< term := (term * x)/n;
result := result + term/n;
n := n + 1>>; >>
else if f = Shi then
<< n:=1; term := x; result := x;
while abs(term) > precis do
<< term := (term * x*x)/(2n * (2n+1));
result := result + term/(2n+1);
n := n + 1>>; >>
else if f = Chi then
<< n:=1; term := 1; result := euler!*constant + log(x);
while abs(term) > precis do
<< term := (term * x*x)/((2n-1) * 2n);
result := result + term/(2n);
n := n + 1>>; >>
else if f = erf then
if x < 0 then result := - compute_int_functions(-x,f) else
<< n:=1; term := x; result := x;
if floor(x*7) > pre then precision floor(x*7);
interm := -1 * x*x;
while abs(term) > precis do
<< term := (term * interm)/n;
result := result + term/(2n+1);
n := n + 1>>;
precision pre;
result := result*2/sqrt(pi) ;>>
else if f = Fresnel_S then
<< if x > 4.0 then precision max(pre,40);
if x > 6.0 then precision max(pre,80);
n:=1; term := x^3*pi/2; result := term/3;
interm := x^4*(pi/2)^2;
while abs(term) > precis do
<< term := -1 * (term * interm)/(2n * (2n+1));
result := result + term/(4n+3);
n := n + 1>>; >>
else if f = Fresnel_C then
<< if x > 4.0 then precision max(pre,40);
if x > 6.0 then precision max(pre,80);
n:=1; term := x; result := x; interm := x^4*(pi/2)^2;
while abs(term) > precis do
<< term := -1 * (term * interm)/(2n * (2n-1));
result := result + term/(4n+1);
n := n + 1>>; >>;
precision pre;
return result
end;
% 16 Jan 96.
algebraic <<
let { binomial (~n,~k) => ((for l:=0:(k-1) product (n-l))/factorial k)
when fixp n and fixp k and k >=0,
binomial (~n,~k) => 1
when fixp n and fixp k and n >= k and k=0,
binomial (~n,~k) => 0
when fixp n and fixp k and n<k and n >=0,
binomial (~n,~k) => 0
when fixp n and fixp k and k < 0,
binomial (~n,~k) => Gamma(n+1) / Gamma (k+1) / Gamma(n-k+1)
when numberp n and numberp k and not(fixp (n - k) and
(n-k) < 0),
df(binomial(~c,~k),c) => binomial(c,k)*(Psi(1+c)-Psi(1+c-k))
} >>;
endpatch;
patch ztrans;
% 27 Nov 95.
algebraic let {binomial(~n,~k)=>
(for i:=0:k-1 product n-i)/factorial(k) when fixp(k)};
endpatch;
endmodule;
end;