module solve; % Solve one or more algebraic equations.
% Author: David R. Stoutemyer.
% Modifications by: Anthony C. Hearn and Donald R. Morrison.
fluid '(!*exp asymplis!*);
global '(!!arbint
!!gcd
!*allbranch
!*micro!-version
!*nonlnr
!*ppsoln
!*solveinterval
!*solvesingular
multiplicities!*);
switch allbranch,solvesingular; % solveinterval.
flag('(multiplicities!*),'share);
% ***** Global Declarations *****
array !!cf(12), !!interval(10,2), !!exact(10);
!*allbranch := t; % Returns all branches of solutions if T;
%!*solveinterval = nil;% Attempts to isolate insoluble, real roots if T;
!*solvesingular := t; % Default value.
% !!gcd SOLVECOEFF returns GCD of powers of its arg in this
% !!cf : Array of coeffs from SOLVECOEFF
algebraic operator arbint, arbreal, intervl, list;
% algebraic operator arbcomplex;
% Done this way since it's also defined in the glmat module.
deflist('((arbcomplex simpiden)),'simpfn);
% ***** Utility Functions *****
symbolic procedure freeofl(u,v);
null v or freeof(u,car v) and freeofl(u,cdr v);
symbolic procedure ratnump x;
% Returns T iff any prefix expression x is a rational number.
atom numr(x := simp!* x) and atom denr x;
flag ('(ratnump), 'direct);
symbolic procedure allkern elst;
% Returns list of all top-level kernels in the list of standard
% forms elst.
if null elst then nil
else union(kernels car numr elst, allkern cdr elst);
symbolic procedure topkern(u,x);
% Returns list of top level kernels in the standard form u that
% contain the kernel x;
for each j in kernels u conc if not freeof(j,x) then list j else nil;
symbolic procedure coeflis ex;
% Ex is a standard form. Returns a list of the coefficients of the
% main variable in ex in the form ((expon . coeff) (expon . coeff)
% ... ), where the expon's occur in increasing order, and entries do
% not occur of zero coefficients.
begin scalar ans,var;
if domainp(ex) then return (0 . ex);
var := mvar(ex);
while (not domainp(ex)) and mvar(ex)=var do
<<ans := (ldeg(ex) . lc(ex)) . ans; ex := red(ex) >>;
if ex then ans := (0 . ex) . ans;
return ans
end;
% ***** Evaluation Interface *****
symbolic procedure solveeval u;
begin scalar arglist; integer nargs;
arglist := u;
nargs := length(arglist);
u := if nargs=1 then solve0(car arglist,nil)
else if nargs=2
then solve0(car arglist, cadr arglist)
else solve0(car arglist,'list . cdr arglist);
return !*solvelist2solveeqlist u
end;
put('solve,'psopfn,'solveeval);
symbolic procedure !*solvelist2solveeqlist u;
begin scalar x,y,z;
for each j in u do
<<if caddr j=0 then rederr "zero multiplicity"
else if null cadr j
then x := for each k in car j collect
list('equal,mk!*sq k,0)
else x := for each k in pair(cadr j,car j)
collect list('equal,car k,mk!*sq cdr k);
if length x > 1 then z := ('list . x) . z
else z := car x . z;
y := caddr j . y>>;
multiplicities!* := 'list . y;
return 'list . z
end;
% ***** Fundamental SOLVE Procedures *****
comment these procedures return the solution of a list of equations as a
list of elements with three fields: the solutions, the variables (or
NIL if the equations could not be solved) and the multiplicity;
symbolic procedure solve0(elst, xlst);
% elst is any prefix expression, including the kernel named LST with
% any number of arguments. XLST is a kernel, perhaps named LIST with
% any number of arguments. Solves eqns in ELST for vars in XLST,
% returning either a list of solutions, or a single solution;
begin scalar !*exp,vars; integer neqn;
!*exp := t;
elst := for each j in solveargchk elst
collect simp!* if eqexpr j then !*eqn2a j else j;
neqn := length elst;
if neqn = 0 then rederr "SOLVE called with no equations";
if null xlst
then <<vars := allkern 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 rederr "SOLVE called with no variables"
else if neqn = 1
then if null numr car elst
then return if !*solvesingular
then list list(list (makearbcomplex() ./ 1),
vars,1)
else nil
else if length vars=1
then return solvesq(car elst,car vars,1);
% more than one equation or variable.
elst := solvesys(for each j in elst collect numr j,vars);
return if null elst then nil
else if null cdr elst then list list(car elst,vars,1)
else if null !*nonlnr then rederr "Unbalanced SOLVE equations"
else elst
end;
symbolic procedure solveargchk u;
if getrtype u eq 'list then cdr reval u
else if atom u or not(car u eq 'lst) then list u
else cdr u;
% ***** Procedures for solving a single eqn *****
symbolic procedure solvesq (ex,var,mul);
% Attempts to find solutions for standard quotient ex with respect to
% top level occurrences of var and kernels containing variable var.
% Solutions containing more than one such kernel are returned
% unsolved, and solve1 is applied to the other solutions. Integer
% mul is the multiplicity passed from any previous factorizations.
% Returns a list of triplets consisting of solutions, variables and
% multiplicity.
begin scalar e1,x1,y,z; integer mu;
ex := numr ex;
if null topkern(ex,var) then return nil;
ex := fctrf ex;
% now process monomial.
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 x1
then z := append(
if null cdr x1 then solve1(e1,car x1,var,mu)
else if (y := principal!-of!-powers!-soln(e1,x1,var,mu))
neq 'unsolved
then y
else if not smemq('sol,
(x1:=simp!* list('sol,mk!*sq(e1 ./ 1), var)))
then solvesq(x1,var,mu)
else list list(list(e1 ./ 1),nil,mu),
z)>>;
return z
end;
symbolic procedure principal!-of!-powers!-soln(ex,x1,var,mu);
% Finds solutions of ex=0 by the principal of powers method, or
% NIL if no such solutions exist.
begin scalar z;
if null !*ppsoln then return 'unsolved;
a: if null x1 then return 'unsolved
else if suitable!-expt car x1
and not((z := pr!-pow!-soln1(ex,car x1,var,mu)) eq 'unsolved)
then return z;
x1 := cdr x1;
go to a
end;
symbolic procedure pr!-pow!-soln1(ex,y,var,mu);
begin scalar oldkord,z;
oldkord := setkorder list y;
z := reorder ex;
setkorder oldkord;
if ldeg z neq 1 then return 'unsolved;
z := coeflis z;
if length z neq 2 or caar z neq 0
then errach list("solve confused",ex,z);
z := exptsq(quotsq(negsq(cdar z ./ 1),cdadr z ./ 1),
caddr caddr y);
z := solvesq(subs2 addsq(simp!* cadr y,negsq z),var,mu);
z := check!-solutions(z,ex);
return z
end;
symbolic procedure check!-solutions(z,ex);
begin scalar x,y;
while z do
if null cadar z then <<z := nil; x := 'unsolved>>
else if null (y := numr subf(ex,list(caadar z .
mk!*sq caaar z)))
or null numvalue y
then <<x := car z . x; z := cdr z>>
else z := cdr z;
return x
end;
symbolic procedure numvalue u;
% Find floating point value of sf u.
begin scalar !*numval,x;
!*numval := t;
x := setdmode('float,t);
u := numr simp prepf u;
if x then setdmode(x,t) else setdmode('float,nil);
return if eqcar(u,'!:ft!:) and 1000000*abs cdr u < 1 then nil
else u
end;
symbolic procedure suitable!-expt u;
eqcar(u,'expt) and eqcar(caddr u,'quotient) and cadr caddr u = 1
and fixp caddr caddr u;
symbolic procedure solve1(e1,x1,var,mu);
comment e1 is a standard form, non-trivial in the kernel x1,
which is itself a function of var, mu is an integer.
Uses roots of unity, known solutions,
inverses, together with quadratic, cubic and quartic
formulas, treating other cases as unsolvable. Returns nil;
begin scalar b,c,coeffs,hipow; integer n;
hipow := errorset(solvecoeff(e1, x1),nil,nil);
if atom hipow then return list list(list(e1 . 1),nil,mu);
% solvecoeff problem - no soln.
hipow := car hipow;
n:= !!gcd; % numerical gcd of powers.
for i := 0:hipow do
coeffs := nilchk getelv list('!!cf,i) . coeffs;
if hipow = 1
then return begin scalar lincoeff,y,z;
b:=prepsq quotsq(negsq cadr coeffs,car coeffs);
if n neq 1 then b := list('expt,b,list('quotient,1,n));
% We may need to merge more solutions in the following if
% there are repeated roots.
for k := 0:n-1 do % equation in power of var.
<<lincoeff := simp!* list('times,b,
mkexp list('quotient,list('times,k,2,'pi),n));
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)) % known 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); % allows for power variable
for each j in apply('solvequadratic,coeffs)
conc solvesq(subtrsq(x1,j),var,mu)>>
else return begin scalar d,f,rcoeffs;
% At this point, we cannot write down the solution directly, so
% we look for various forms that we know how to solve.
f:=(hipow+1)/2;
d:=exptsq(simp!* x1,n);
rcoeffs := reverse coeffs;
return if solve1test1(coeffs,rcoeffs,f) % coefficients symmetric
then if f+f=hipow+1 % odd
then <<c:=addsq(d, 1 ./ 1);
append(solvesq(c,var,mu),
solvesq(quotsq(e1 ./ 1, c),var,mu))>>
else <<setelv(list('!!cf,0),2 ./ 1);
setelv(list('!!cf, 1), simp!* '!!x);
c:=addsq(multsq(getelv(list('!!cf,f+1)),
getelv(list('!!cf,1))),
getelv(list('!!cf,f)));
for j:=2:f do <<
setelv(list('!!cf, j),
subtrsq(multsq(getelv(list('!!cf,1)),
getelv(list('!!cf,j-1))),
getelv(list('!!cf,j-2))));
c:=addsq(c,multsq(getelv(list('!!cf,j)),
getelv(list('!!cf,f+j))))>>;
for each j in solvesq(c,'!!x,mu) conc
solvesq(addsq(1 ./ 1,multsq(d,subtrsq(d,caar j))),
var,caddr j)>>
else if solve1test2(coeffs,rcoeffs,f)
% coefficients antisymmetric
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)>>
% equation has no symmetry
else if hipow=3 and null !*micro!-version
then for each j in apply('solvecubic,coeffs)
conc solvesq(subtrsq(d,j),var,mu)
else if hipow=4 and null !*micro!-version
then for each j in apply('solvequartic,coeffs)
conc solvesq(subtrsq(d,j),var,mu)
else if !*solveinterval and univariatep e1
then solveinterval(e1,var,mu)
else list list(list(e1 ./ 1),nil,mu)
% We can't solve quintic and higher
end
end;
symbolic procedure solnmerge(u,varlist,mu,y);
% Merge solutions in case of multiplicities. It may be that this is
% only needed for the trivial solution x=0.
if null y then list list(u,varlist,mu)
else if u = caar y and varlist = cadar y
then list(caar y,cadar y,mu+caddar y) . cdr y
else car y . solnmerge(u,varlist,mu,cdr y);
symbolic procedure nilchk u; if null u then !*f2q u else u;
symbolic procedure solve1test1(coeffs,rcoeffs,f);
% True if equation is symmetric in its coefficients. f is midpoint.
begin integer j;
a: if j>f then return t
else if car coeffs neq car rcoeffs then return nil;
coeffs := cdr coeffs;
rcoeffs := cdr rcoeffs;
j := j+1;
go to a
end;
symbolic procedure solve1test2(coeffs,rcoeffs,f);
% True if equation is antisymmetric in its coefficients. f is
% midpoint.
begin integer j;
a: if j>f then return t
else if numr addsq(car coeffs,car rcoeffs) then return nil;
coeffs := cdr coeffs;
rcoeffs := cdr rcoeffs;
j := j+1;
go to a
end;
symbolic procedure solveabs u;
begin scalar mu,var,lincoeff;
var := cadr u;
mu := caddr u;
lincoeff := cadddr u;
u := simp!* caar u;
return append(solvesq(addsq(u,lincoeff),var,mu),
solvesq(subtrsq(u,lincoeff),var,mu))
end;
put('abs,'solvefn,'solveabs);
symbolic procedure solveexpt u;
begin scalar c,mu,var,lincoeff;
var := cadr u;
mu := caddr u;
lincoeff := cadddr u;
u := car u;
return if freeof(car u,var) % c**(...) = b.
then <<if !*allbranch
then <<!!arbint:=!!arbint+1;
c:=list('times,2,'i,'pi,
list('arbint,!!arbint))>>
else c:=0;
solvesq(subtrsq(simp!* cadr u,
quotsq(addsq(simp!* list('log,mk!*sq lincoeff),
simp!* c),
simp!* list('log,car u))),var,mu)>>
else if freeof(cadr u,var) % (...)**(m/n) = b;
then if ratnump cadr u
then solve!-fractional!-power(u,lincoeff,var,mu)
else << % (...)**c = b.
if !*allbranch
then <<!!arbint:=!!arbint+1;
c := mkexp list('times,
list('arbreal,!!arbint))>>
else c:=1;
solvesq(subtrsq(simp!* car u,
multsq(simp!* list('expt,
mk!*sq lincoeff,
mk!*sq invsq
simp!* cadr u),
simp!* c)),var,mu)>>
% (...)**(...) = b : transcendental.
else list list(list subtrsq(simp!*('expt . u),lincoeff),nil,mu)
end;
symbolic procedure solve!-fractional!-power(u,x,var,mu);
% attempts solution of equation car u**cadr u=x with respect to
% kernel var and with multiplicity mu, where cadr u is a rational
% number.
begin scalar v,w,z;
v := simp!* car u;
w := simp!* cadr u;
z := solvesq(subs2 subtrsq(exptsq(v,numr w),exptsq(x,denr w)),
var,mu);
w := subtrsq(simp('expt . u),x);
z := check!-solutions(z,numr w);
return if z eq 'unsolved then list list(list w,nil,mu) else z
end;
put('expt,'solvefn,'solveexpt);
symbolic procedure solvelog u;
solvesq(subtrsq(simp!* caar u,simp!* list('expt,'e,mk!*sq cadddr u)),
cadr u,caddr u);
put('log,'solvefn,'solvelog);
symbolic procedure solvecos u;
begin scalar c,d,z;
if !*allbranch
then <<!!arbint:=!!arbint+1;
c:=list('times,2,'pi,list('arbint,!!arbint))>>
else c:=0;
c:=subtrsq(simp!* caar u,simp!* c);
d:=simp!* list('acos,mk!*sq cadddr u);
z := solvesq(subtrsq(c,d), cadr u,caddr u);
if !*allbranch
then z := append(solvesq(addsq(c,d), cadr u,caddr u),z);
return z
end;
put('cos,'solvefn,'solvecos);
symbolic procedure solvesin u;
begin scalar c,d,f,z;
if !*allbranch
then <<!!arbint:=!!arbint+1;
f:=list('times,2,'pi,list('arbint,!!arbint))>>
else f:=0;
c:=simp!* caar u;
d:=list('asin,mk!*sq cadddr u);
z := solvesq(subtrsq(c,simp!* list('plus,d,f)),cadr u,caddr u);
if !*allbranch
then z := append(solvesq(subtrsq(c,simp!* list('plus,'pi,
mk!*sq subtrsq(simp!* f,simp!* d))),
cadr u,caddr u),z);
return z
end;
put('sin,'solvefn,'solvesin);
symbolic procedure mkexp u;
list('plus,list('cos,x),list('times,'i,list('sin,x)))
where x = reval u;
symbolic procedure solvecoeff(ex,var);
% ex is a standard form and var a kernel. Puts the coefficients
% (as standard quotients) of var in ex into the elements of !!cf,
% with index equal to the exponent divided by the gcd of all the
% exponents. This GCD is put into !!GCD, and the highest power
% divided by the gcd is put into hipow. Returns hipow. Note that
% !!cf (an array), !!gcd a global.
begin scalar clist,hipow,oldkord;
oldkord := setkorder list var;
clist := reorder ex;
setkorder oldkord;
hipow := ldeg clist;
clist := coeflis clist;
!!gcd := caar clist;
for each x in cdr clist do !!gcd := gcdn(car x, !!gcd);
for i := 0:(car(dimension('!!cf))-1) do setelv(list('!!cf, i), nil);
for each x in clist do setelv(list('!!cf, car x/!!gcd),cdr x ./ 1);
hipow := hipow/!!gcd;
return hipow
end;
symbolic procedure solveinterval(ex,var,mu);
% ex is a standard form, var the relevant variable and mu the root
% multiplicity. Isolates insoluble, real roots of EX in rational
% intervals, returning solutions in terms of INTERVL(Lowlim,Highlim).
begin scalar z;
realroot(prepf ex,prepsq !*k2q mvar ex,'!!interval,'!!exact);
for i := 1:getelv list('!!exact,0) do
z := list(list simp!* getelv list('!!exact,i),list var,mu) . z;
for i := 1:getelv list('!!interval,0,0) do
z := list(list simp!* list('intervl,
getelv list('!!interval,i,1),
getelv list('!!interval,i,2)),
list var,mu). z;
return z
end;
symbolic procedure realroot(u,v,w,x);
rederr("Real root finding not yet implemented");
% ***** Procedures for solving a system of eqns *****
symbolic procedure solvesys(exlist,varlis);
% exlist is a list of standard forms, varlis a list of kernels. If
% the elements of varlis are linear in the elements of exlist, and
% further the system of linear eqns so defined is non-singular, then
% SOLVESYS returns a list of a list of standard quotients which are
% solutions of the system, ordered as in varlis. Otherwise an error
% results.
begin scalar eqtype,oldkord;
oldkord := setkorder varlis;
exlist := for each j in exlist collect reorder j;
% See if equations are linear or non-linear.
eqtype := 'solvelnrsys;
for each ex in exlist do
for each var in varlis do
if not domainp ex and mvar ex=var
then if ldeg ex>1 or not freeofl(lc ex,varlis)
then eqtype := 'solvenonlnrsys
else ex := red ex;
if eqtype eq 'solvenonlnrsys and null !*nonlnr
then rederr "Non linear equation solving not yet implemented";
exlist:=errorset(list(eqtype,mkquote exlist,mkquote varlis),t,t);
setkorder oldkord;
if errorp exlist then error1() else return car exlist
end;
endmodule;
module glsolve; % Routines for solving a general system of linear eqns.
% Author: Eberhard Schruefer.
%**********************************************************************
%*** The number of equations and the number of unknowns are ***
%*** arbitrary i.e. the system can be under- or overdetermined. ***
%*** Method used is Cramer's rule, realized through exterior ***
%*** multiplication. ***
%**********************************************************************
fluid '(kord!*);
global '(!!arbint !*solvesingular);
% algebraic operator arbcomplex; % Already defined in main solve module.
symbolic procedure glsolve!-eval(u,bool);
% This allows glsolve to be called at the user level. I'm not
% sure this is now a good idea, since this code does not check
% for non-linear equations and so on.
begin scalar unknowns,equations,okord,solutions;
if cdr u then
unknowns := for each j in cdadr u collect !*a2k j;
okord := setkorder append(unknowns,kord!*);
equations := for each j in cdar u collect
reorder numr simp!* j;
if null unknowns then unknowns := allkernf equations;
solutions := glnrsolve(equations,unknowns);
setkorder okord;
if null solutions then return '(list); % empty list.
solutions := nil . solutions;
return 'list .
for each j in unknowns collect
list('equal,j,mk!*sq car(solutions := cdr solutions))
end;
symbolic procedure allkernf u;
if null u then nil else union(kernels car u,allkernf cdr u);
put('glsolve,'psopfn,'glsolve!-eval);
symbolic procedure solvelnrsys(u,v);
% This is hook to general solve package. u is a list of polynomials
% (s.f.'s) linear in the kernels of list v. Result is a matrix
% standard form for the solutions.
list glnrsolve(u,v);
symbolic procedure glnrsolve(u,v);
%u is a list of polynomials (s.f.'s) linear in the kernels
%of list v. Result is an untagged list of solutions.
begin scalar arbvars,sgn,x,y;
x := !*sf2ex(car u,v);
u := cdr u;
for each j in u do
if y := extmult(!*sf2ex(j,v),x)
then x := y;
if inconsistency!-chk x
then rederr "SOLVE given inconsistent equations";
arbvars := for each j in setdiff(v,lpow x) collect
j . makearbcomplex();
if arbvars and null !*solvesingular
then rederr "SOLVE given singular equations";
if null red x then return
for each j in v collect
if y := atsoc(j,arbvars) then !*f2q cdr y else nil ./ 1;
sgn := evenp length lpow x;
return for each j in v collect if y := atsoc(j,arbvars)
then !*f2q cdr y
else mkglsol(j,x,sgn := not sgn,arbvars)
end;
symbolic procedure inconsistency!-chk u;
null u or ((nil memq lpow u) and inconsistency!-chk red u);
symbolic procedure mkglsol(u,v,sgn,arbvars);
begin scalar s,x,y;
x := nil ./ 1;
y := lpow v;
for each j on red v do
if s := glsolterm(u,y,j,arbvars)
then x := addsq(cancel(s ./ lc v),x);
return if sgn then negsq x else x
end;
symbolic procedure glsolterm(u,v,w,arbvars);
begin scalar x,y,sgn;
x := lpow w;
a: if null x then return
if null car y then lc w
else multf(cdr atsoc(car y,arbvars),
if sgn then negf lc w else lc w);
if car x eq u then return nil
else if car x memq v then <<x := cdr x;
if y then sgn := not sgn>>
else if y then return nil
else <<y := list car x; x := cdr x>>;
go to a
end;
%**** Support for exterior multiplication ****
% Data structure is lpow ::= list of variables in exterior product
% lc ::= standard form
symbolic procedure !*sf2ex(u,v);
%Converts standardform u into a form distributed w.r.t. v
%*** Should we check here if lc is free of v?
if null u then nil
else if domainp u or null(mvar u memq v) then list nil .* u .+ nil
else list mvar u .* lc u .+ !*sf2ex(red u,v);
symbolic procedure extmult(u,v);
%Special exterior multiplication routine. Degree of form v is
%arbitrary, u is a one-form.
if null u or null v then nil
else (if x then cdr x .* (if car x then negf multf(lc u,lc v)
else multf(lc u,lc v))
.+ extadd(extmult(!*t2f lt u,red v),
extmult(red u,v))
else extadd(extmult(red u,v),extmult(!*t2f lt u,red v)))
where x = ordexn(car lpow u,lpow v);
symbolic procedure extadd(u,v);
if null u then v
else if null v then u
else if lpow u = lpow v then
(lambda x,y; if null x then y else lpow u .* x .+ y)
(addf(lc u,lc v),extadd(red u,red v))
else if ordexp(lpow u,lpow v) then lt u .+ extadd(red u,v)
else lt v .+ extadd(u,red v);
symbolic procedure ordexp(u,v);
if null u then t
else if car u eq car v then ordexp(cdr u,cdr v)
else if null car u then nil
else if null car v then t
else ordop(car u,car v);
symbolic procedure ordexn(u,v);
%u is a single variable, v a list. Returns nil if u is a member
%of v or a dotted pair of a permutation indicator and the ordered
%list of u merged into v.
begin scalar s,x;
a: if null v then return(s . reverse(u . x))
else if u eq car v then return nil
else if u and ordop(u,car v) then
return(s . append(reverse(u . x),v))
else <<x := car v . x;
v := cdr v;
s := not s>>;
go to a
end;
endmodule;
module quartic; % Procedures for solving cubic, quadratic and quartic
% eqns.
% Author: Anthony C. Hearn.
fluid '(!*sub2);
symbolic procedure multfq(u,v);
% Multiplies standard form U by standard quotient V.
begin scalar x;
x := gcdf(u,denr v);
return multf(quotf(u,x),numr v) ./ quotf(denr v,x)
end;
symbolic procedure quotsqf(u,v);
% Forms quotient of standard quotient U and standard form V.
begin scalar x;
x := gcdf(numr u,v);
return quotf(numr u,x) ./ multf(quotf(v,x),denr u)
end;
symbolic procedure cubertq u;
simpexpt list(mk!*sq subs2!* u,'(quotient 1 3));
% SIMPRAD(U,3);
symbolic procedure sqrtq u;
simpexpt list(mk!*sq subs2!* u,'(quotient 1 2));
% SIMPRAD(U,2);
symbolic procedure subs2!* u; <<!*sub2 := t; subs2 u>>;
symbolic procedure solvequadratic(a2,a1,a0);
% a2, a1 and a0 are standard quotients.
% solves a2*x**2+a1*x+a0=0 for x.
% returns a list of standard quotient solutions.
begin scalar d;
d := sqrtq subtrsq(quotsqf(exptsq(a1,2),4),multsq(a2,a0));
a1 := quotsqf(negsq a1,2);
return list(subs2!* quotsq(addsq(a1,d),a2),
subs2!* quotsq(subtrsq(a1,d),a2))
end;
symbolic procedure solvecubic(a3,a2,a1,a0);
% a3, a2, a1 and a0 are standard quotients.
% solves a3*x**3+a2*x**2+a1*x+a0=0 for x.
% returns a list of standard quotient solutions.
% See Abramowitz and Stegun, Sect. 3.8.2, for details.
begin scalar q,r,sm,sp,s1,s2,x;
a2 := quotsq(a2,a3);
a1 := quotsq(a1,a3);
a0 := quotsq(a0,a3);
q := subtrsq(quotsqf(a1,3),quotsqf(exptsq(a2,2),9));
r := subtrsq(quotsqf(subtrsq(multsq(a1,a2),multfq(3,a0)),6),
quotsqf(exptsq(a2,3),27));
x := sqrtq addsq(exptsq(q,3),exptsq(r,2));
s1 := cubertq addsq(r,x);
s2 := if numr s1 then negsq quotsq(q,s1)
else cubertq subtrsq(r,x);
% This optimization only works if s1 is non zero.
sp := addsq(s1,s2);
sm := quotsqf(multsq(simp '(times i (sqrt 3)),subtrsq(s1,s2)),2);
x := subtrsq(sp,quotsqf(a2,3));
sp := negsq addsq(quotsqf(sp,2),quotsqf(a2,3));
return list(subs2!* x,subs2!* addsq(sp,sm),
subs2!* subtrsq(sp,sm))
end;
symbolic procedure solvequartic(a4,a3,a2,a1,a0);
% Solve the quartic equation a4*x**4+a3*x**3+a2*x**2+a1*x+a0 = 0,
% where the ai are standard quotients, using technique described in
% Section 3.8.3 of Abramowitz and Stegun;
begin scalar x,y,z;
% Convert equation to monomial form.
a3 := quotsq(a3,a4);
a2 := quotsq(a2,a4);
a1 := quotsq(a1,a4);
a0 := quotsq(a0,a4);
% Build and solve the resultant cubic equation. We select an
% arbitrary member of its set of solutions. Ideally we should
% only generate one solution, which should be the simplest.
y := subtrsq(exptsq(a3,2),multfq(4,a2));
% note that only first cubic solution is used here. We could save
% computation by using this fact.
x := car solvecubic(!*f2q 1,
negsq a2,
subs2!* subtrsq(multsq(a1,a3),multfq(4,a0)),
subs2!* negsq addsq(exptsq(a1,2),
multsq(a0,y)));
% Now solve the two equivalent quadratic equations.
y := sqrtq addsq(quotsqf(y,4),x);
z := sqrtq subtrsq(quotsqf(exptsq(x,2),4),a0);
a3 := quotsqf(a3,2);
x := quotsqf(x,2);
return append(solvequadratic(!*f2q 1,addsq(a3,y),subtrsq(x,z)),
solvequadratic(!*f2q 1,subtrsq(a3,y),addsq(x,z)))
end;
endmodule;
module solvetab; % Simplification rules for SOLVE.
% Author: David R. Stoutemyer.
% Modifications by: Anthony C. Hearn and Donald R. Morrison;
algebraic operator sol;
put('asin, 'inverse, 'sin);
put('acos, 'inverse, 'cos);
algebraic;
comment Rules for reducing the number of distinct kernels in an
equation;
for all a,b,x such that ratnump c and ratnump d let
sol(a**c-b**d, x) = a**(c*lcm(c,d)) - b**(d*lcm(c,d));
for all a,b,c,d,x such that a freeof x and c freeof x let
sol(a**b-c**d, x) = e**(b*log a - d*log c);
for all a,b,c,d,x such that a freeof x and c freeof x let
sol(a*log b + c*log d, x) = b**a*d**c - 1,
sol(a*log b - c*log d, x) = b**a - d**c;
for all a,b,c,d,f,x such that a freeof x and c freeof x let
sol(a*log b + c*log d + f, x) = sol(log(b**a*d**c) + f, x),
sol(a*log b + c*log d - f, x) = sol(log(b**a*d**c) - f, x),
sol(a*log b - c*log d + f, x) = sol(log(b**a/d**c) + f, x),
sol(a*log b - c*log d - f, x) = sol(log(b**a/d**c) - f, x);
for all a,b,d,f,x such that a freeof x let
sol(a*log b + log d + f, x) = sol(log(b**a*d) + f, x),
sol(a*log b + log d - f, x) = sol(log(b**a*d) - f, x),
sol(a*log b - log d + f, x) = sol(log(b**a/d) + f, x),
sol(a*log b - log d - f, x) = sol(log(b**a/d) - f, x),
sol(log d - a*log b + f, x) = sol(log(d/b**a) + f, x),
sol(log d - a*log b - f, x) = sol(log(d/b**a) - f, x);
for all a,b,c,d,x such that a freeof x and c freeof x let
sol(a*log b + c*log d, x) = b**a*d**c - 1,
sol(a*log b - c*log d, x) = b**a - d**c;
for all a,b,d,x such that a freeof x let
sol(a*log b + log d, x) = b**a*d - 1,
sol(a*log b - log d, x) = b**a - d,
sol(log d - a*log b, x) = d - b**a;
for all a,b,c,x let
sol(log a + log b + c, x) = sol(log(a*b) + c, x),
sol(log a - log b + c, x) = sol(log(a/b) + c, x),
sol(log a + log b - c, x) = sol(log(a*b) - c, x),
sol(log a - log b - c, x) = sol(log(a/b) - c, x);
for all a,c,x such that c freeof x let
sol(log a + c, x) = a - e**c,
sol(log a - c, x) = a - e**(-c);
for all a,b,x let
sol(log a + log b, x) = a*b - 1,
sol(log a - log b, x) = a - b,
sol(cos a - sin b, x) = sol(cos a - cos(pi/2-b), x),
sol(sin a + cos b, x) = sol(sin a - sin(b-pi/2), x),
sol(sin a - cos b, x) = sol(sin a - sin(pi/2-b), x),
sol(sin a + sin b, x) = sol(sin a - sin(-b), x),
sol(sin a - sin b, x) = if !*allbranch then sin((a-b)/2)*
cos((a+b)/2) else a-b,
sol(cos a + cos b, x) = if !*allbranch then cos((a+b)/2)*
cos((a-b)/2) else a+b,
sol(cos a - cos b, x) = if !*allbranch then sin((a+b)/2)*
sin((a-b)/2) else a-b,
sol(asin a - asin b, x) = a-b,
sol(asin a + asin b, x) = a+b,
sol(acos a - acos b, x) = a-b,
sol(acos a + acos b, x) = a+b;
symbolic;
endmodule;
end;