module solvealg; % Solution of equations and systems which can
% be lifted to algebraic (polynomial) systems.
% Author: Herbert Melenk.
% Copyright (c) 1992 The RAND Corporation and Konrad-Zuse-Zentrum.
% All rights reserved.
% August 1992: added material for
% rule set for reduction of trig. polynomial terms to
% elementary expressions in sin and cos,
% constant expressions in sin, cos and constant roots,
% closed form results for trigonometric systems.
% general exponentials.
% avoiding false solutions with surds.
%
% May 1993: better handling of products of exponentials
% with common base,
% additional computation branch for linear parts of
% nonlinear systems.
% July 1996: safe handling of twice (or more) the same input
% (not handling the case, that one equation is a multiple
% of an other one)
fluid '(!*expandexpt); % from simp.red
fluid '( system!* % system to be solved
osystem!* % original system on input
uv!* % user supplied variables
iv!* % internal variables
fv!* % restricted variables
kl!* % kernels to be investigated
sub!* % global substitutions
inv!* % global inverse substitutions
depl!* % REDUCE dependency list
!*solvealgp % true if using this module
solvealgdb!* % collecting some data
last!-vars!* % collection of innermost aux variables
const!-vars!* % variables representing constants
root!-vars!* % variables representing root expressions
!*expli % local switch: explicit solution
groebroots!* % predefined roots from input surds
!*test_solvealg % debugging support
!*arbvars
!*varopt
solve!-gensymcounter
);
fluid '(!*trnonlnr);
% If set on, the modified system and the Groebner result
% or the reason for the failure are printed.
global '(loaded!-packages!* !!arbint);
switch trnonlnr;
!*solvealgp := t;
% Solvenonlnrsys receives a system of standard forms and
% a list of variables from SOLVE. The system is lifted to
% a polynomial system (if possible) in substituting the
% non-atomic kernels by new variables and appending additonal
% relations, e.g.
% replace add
% sin u,cos u -> su,cu su^2+cu^2-1
% u^(1/3) -> v v^3 - u
% ...
% in a recursive style. If completely successful, the
% system definitely can be treated by Groebner or any
% other polynomial system solver.
%
% Return value is a pair
% (tag . res)
% where "res" is nil or a structure for !*solvelist2solveeqlist
% and "tag" is one of the following:
%
% T a satisfactory solution was generated,
%
% FAILED the algorithm cannot be applied (res=nil)
%
% INCONSISTENT the algorithm could prove that the
% the system has no solution (res=nil)
%
% NIL the complexity of the system could
% be reduced, but some (or all) relations
% remain still implicit.
% rules to be applied locally for converting composite transcendental
% function forms into simpler ones
algebraic <<
solvealg!-rules1:=
{ sin(~alpha + ~beta) => sin(alpha)*cos(beta) + cos(alpha)*sin(beta),
cos(~alpha + ~beta) => cos(alpha)*cos(beta) - sin(alpha)*sin(beta),
sin(~n*~alpha) => sin(alpha)*cos((n-1)*alpha)
+ cos(alpha)*sin((n-1)*alpha) when fixp n,
cos(~n*~alpha) => cos(alpha)*cos((n-1)*alpha)
- sin(alpha)*sin((n-1)*alpha) when fixp n,
sin(~alpha)**2 => 1 - cos(alpha)**2,
sinh(~alpha+~beta) => sinh(alpha)*cosh(beta) + cosh(alpha)*sinh(beta),
cosh(~alpha+~beta) => cosh(alpha)*cosh(beta) + sinh(alpha)*sinh(beta),
sinh(~n*~alpha) => sinh(alpha)*cosh((n-1)*alpha)
+ cosh(alpha)*sinh((n-1)*alpha) when fixp n,
cosh(~n*~alpha) => cosh(alpha)*cosh((n-1)*alpha)
+ sinh(alpha)*sinh((n-1)*alpha) when fixp n,
sinh(~alpha)**2 => cosh(alpha)**2 - 1};
solvealg!-rules2:=
{ tan(~alpha) => sin(alpha)/cos(alpha),
cot(~alpha) => cos(alpha)/sin(alpha),
tanh(~alpha) => sinh(alpha)/cosh(alpha),
coth(~alpha) => cosh(alpha)/sinh(alpha) } ;
solvealg!-rules3:=
{ sin(~alpha)**2 => 1 - cos(alpha)**2,
sinh(~alpha)**2 => cosh(alpha)**2 - 1};
% Artificial operator for matching powers in a
% product.
operator my!-expt;
solvealg!-rules4:=
{ my!-expt(~a,~b)*my!-expt(a,~c) => my!-expt(a,b+c),
my!-expt(~a,~b)*a => my!-expt(a,b+1)
% my!-expt(~a,~b)/my!-expt(a,~c) => my!-expt(a,b-c)
};
>>;
symbolic procedure solvenonlnrsys(sys,uv);
% interface to algebraic system solver.
% factorize the system and collect solutions.
% After factoring we resimplify with *expandexpt off
% in order to have exponentials to one basis
% collected.
begin scalar q,r,s,tag,!*expandexpt;
s := sys; sys := nil;
for each x in s do sys := union(sys,{x});
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
<<% collect exponentials with same base.
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));
symbolic procedure solvenonlnrcollectexpt u;
% u is a list of standard forms. Reform these
% such that products of exponentials with same basis
% are collected.
if not smemq('expt,u) then u else
<<eval'(let0 '(solvealg!-rules4));
u:=for each q in u collect
numr simp subst('expt,'my!-expt,
reval prepf subst('my!-expt,'expt,q));
eval'(clearrules '(solvealg!-rules4));
u>>;
symbolic procedure solvenonlnrsyslin(eqs,vars,mode);
% Eqs is a system of equations (standard forms,
% implicitly equated to zero); this routine tries
% to reduce the system recursively by separation,
% if one variable occurs in one equation only linearly.
% Mode=NIL: simple version: only pure linear variables
% are substituted.
% T: extended version: replacing variables with
% degree 1 and potentially complicated
% coefficients.
% Returns solution or
% '(not) if not applicable
% '(failed) if applicable but solution failed.
begin scalar d,e,e1,n,s,q,x,v,w,w1,neqs,nvars;
v:=vars;
var_loop:
if null v then return '(not);
x:=car v; v:=cdr v;
w:=eqs;
eqn_loop:
if null w then goto var_loop;
e:=car w; w:=cdr w;
if null e then goto eqn_loop;
if domainp e then return '(inconsistent);
e1:= reorder e where kord!*={x};
if not(mvar e1 =x) or ldeg e1>1 or
smemq(x,d:=lc e1) or smemq(x,n:=red e1) then goto eqn_loop;
if not mode then
<<w:=nil;
for each y in vars do w:=w or smemq(y,d);
if w then return '(not); >>;
% linear form found: n*x+d=0. This is basis for a solution
% x=-n/d. In a second branch the case {n=0,d=0} has to
% be considered if n and d are not constants.
n := reorder n; d:=reorder d;
% Step 1: substitute in remaining equations, solve
% and add linear formula to result.
s:= quotsq(negf n ./ 1, d ./ 1);
neqs := for each eqn in delete(e,eqs) join
<<q:=numr subf(eqn,{x.prepsq s}); if q then {q}>>;
nvars:=for each y in delete(x,vars) join
if smemq(y,neqs) then {y};
w:= if null neqs then '(t (nil nil 1)) else
if null nvars then '(inconsistent) else
if cdr neqs then solvenonlnrsys(neqs,nvars)
else solvenonlnrsysone(car neqs,car nvars);
if car w eq 'failed then return w;
w:=add!-variable!-to!-tagged!-solutions(x,s,w);
% Step 2: add an eventual solution for n=0,d=0;
if domainp d or not mode then return w;
w1:=solvenonlnrsys(n.d.eqs,vars);
return merge!-two!-tagged!-solutions(w,w1);
end;
symbolic procedure solvenonlnrsysone(f,x);
% equation system has been reduced to one. Using solvesq.
begin scalar w;
w:=solvesq(f ./ 1,x,1);
if null w then return '(inconsistent)
else if null cadr car w then return '(failed);
% if not smemq('root_of,w) then goto ret;
% % here we try to find out whether a root_of
% % expression is a useful information or whether
% % it is simply an echo of the input.
% if cdr w then goto ret; % multiple branches: good.
% q := prepsq caar car w;
% if not eqcar(q,'root_of) % not on top level: good.
% then goto ret;
% q:=subst(x,caddr q,cadr q);
% if f = numr simp q then return '(failed);
%ret:
return t.w;
end;
symbolic procedure add!-variable!-to!-tagged!-solutions(x,s,y);
% Y is a tagged solution. Add equation x=s to all members.
if eqcar(y,'inconsistent) then y else
if null y or null cdr y then {t,{{s},{x},1}} else
car y . for each q in cdr y collect
% Put new solution into the last position.
{append(car q,{s}),append(cadr q,{x}),caddr q};
symbolic procedure merge!-two!-tagged!-solutions(w1,w2);
% w1 and w2 are tagged solution sets. Merge these and
% eliminated inconsistent cases.
if car w1='failed or car w2='failed then '(failed) else
if car w1='inconsistent then w2 else
if car w2='inconsistent then w1 else
car w1 . append(cdr w1,cdr w2);
symbolic procedure solvenonlnrsyssep(eqs,vars);
% Eqs is a system of equations (standard forms,
% implicitly equated to zero); this routine tries
% to reduce the system recursively by separation,
% if one variable occurs only in one equation.
begin scalar y,r,s,r0,u,w,tag;
if null vars then return '(failed) else
if null cdr eqs then
<<if not smember(car vars,car eqs) then
return solvenonlnrsyssep(eqs,cdr vars);
r:=solvesq(!*f2q car eqs,car vars,1);
return if r and cadr car r then 't.r else '(failed);
>>;
for each x in vars do if null y then
<<r:=nil;
for each u in eqs do if smember(x,u) then
r:=u.r;
if r and null cdr r then y:=x;
>>;
if null y then return '(failed);
r:=car r;
s:=solvenonlnrsys(delete(r,eqs),delete(y,vars));
if car s='failed then return s else s:=cdr s;
tag := t;
u:=for each s0 in s join
<<
w:=for each q in pair(cadr s0,car s0) join
if not smemq('root_of,cdr q) then
{car q.prepsq cdr q};
r0:=subf(r,w);
r0:=solvesq(r0,y,caddr s0);
if null r0 or null cadr car r0 then tag:='failed;
for each r1 in r0 collect
{caar r1. car s0,y.cadr s0,caddr r1}
>>;
return tag.u;
end;
symbolic procedure solve!-psysp(s,uv);
% T if s is a pure polynomial system.
null s or (solve!-psysp1(car s,uv) and solve!-psysp(cdr s,uv));
symbolic procedure solve!-psysp1(f,uv);
domainp f or
((member(mvar f,uv) or solve!-psysp2(mvar f,uv))
and solve!-psysp1(lc f,uv) and solve!-psysp1(red f,uv));
symbolic procedure solve!-psysp2(v,uv);
% t if there is no interaction between v and uv.
null uv or (not smember(car uv,v) and solve!-psysp2(v,cdr uv));
symbolic procedure solvenonlnrsys1(system!*,uv!*);
% solve one system.
begin scalar r,rules;
osystem!* := system!*;
if solvealgtrig0 system!* then rules:='(solvealg!-rules1);
if smemq('tan,system!*) or smemq('cot,system!*) or
smemq('tanh,system!*) or smemq('coth,system!*) then
rules:='solvealg!-rules2.rules;
r := evalletsub2({rules,'(solvenonlnrsyspre)},nil);
if errorp r then return '(failed) else system!* := car r;
r := solvenonlnrsys2();
return r;
end;
symbolic procedure solvenonlnrsyspre();
(for each p in system!* collect numr simp prepf p)
where dmode!* = nil;
symbolic procedure solvenlnrsimp(u);
% a prepsq including resimplification with additional rules.
% begin scalar r;
% r := evalletsub2({'(solvealg!-rules3),
% {'simp!* ,mkquote u}},nil);
% if errorp r then error(99,list("error during postprocessing simp"));
% return car r;
% end;
simp!* u;
symbolic procedure solvenonlnrsys2();
% Main driver. We need non-local exits here
% because of possibly hidden non algebraic variable
% dependencies.
if null !*solvealgp then system!*:='(FAILED) else % against recursion.
(begin scalar iv!*,kl!*,inv!*,fv!*,r,w,!*solvealgp,solvealgdb!*,sub!*;
scalar last!-vars!*,groebroots!*,const!-vars!*,root!-vars!*;
% preserving the variable sequence if *varopt is off
if not !*varopt then depl!* :=
append(pair(uv!*,append(cdr uv!*,{gensym()})),depl!*);
% hiding dmode because exponentials need integers.
for each f in system!* do solvealgk0
(if dmode!* then numr subf(f,nil) where dmode!*=nil else f);
if !*trnonlnr then print list("original kernels:",kl!*);
if null cdr system!* then
if (smemq('sin,system!*)or smemq('cos,system!*)) and
(r:=solvenonlnrtansub(prepf(w:=car system!*),car uv!*))
and car r
then return solvenonlnrtansolve(r,car uv!*,w)
else if (smemq('sinh,system!*)or smemq('cosh,system!*)) and
(r:=solvenonlnrtanhsub(prepf(w:=car system!*),car uv!*))
and car r
then return solvenonlnrtanhsolve(r,car uv!*,w);
if atom (errorset('(solvealgk1),!*trnonlnr,nil)) where dmode!*=nil
then return (system!*:='(FAILED));
system!*:='LIST.for each p in system!* collect prepf p;
if not('groebner memq loaded!-packages!*)
then load!-package 'groebner;
for each x in iv!* do if not member(x,last!-vars!*) then
for each y in last!-vars!* do depend1(x,y,t);
iv!* := sort(iv!*,function (lambda(a,b);depends(a,b)));
if !*trnonlnr then
<< prin2t "Entering Groebner for system";
writepri(mkquote system!*,'only);
writepri(mkquote('LIST.iv!*), 'only);
>>;
r := list(system!*,'LIST.iv!*);
r := groesolveeval r;
if !*trnonlnr then
<< prin2t "leaving Groebner with intermediate result";
writepri(mkquote r,'only);
terpri(); terpri();
>>;
if 'sin memq solvealgdb!* then r:=solvealgtrig2 r;
if 'sinh memq solvealgdb!* then r:=solvealghyp2 r;
r:= if r='(LIST) then '(INCONSISTENT) else solvealginv r;
system!* := r; % set value aside
return r;
end) where depl!*=depl!* ;
symbolic procedure solvealgk0(p);
% Extract new top level kernels from form p.
if domainp p then nil else
<<if not member(mvar p,kl!*)
and not member(mvar p,iv!*)
then kl!*:=mvar p.kl!*;
solvealgk0(lc p);
solvealgk0(red p)>>;
symbolic procedure solvealgk1();
% Process all kernels in kl!*. Note that kl!* might
% change during processing.
begin scalar k,kl0,kl1;
k := car kl!*;
while k do
<<kl0 := k.kl0;
solvealgk2(k);
kl1 := kl!*; k:= nil;
while kl1 and null k do
if not member(car kl1,kl0) then k:=car kl1
else kl1:=cdr kl1;
>>;
end;
symbolic procedure solvealgk2(k);
% Process one kernel.
(if member(k,uv!*) then solvealgvb0 k and (iv!*:= k.iv!*) else
if atom k then t else
if eq(car k,'EXPT) then solvealgexpt(k,x) else
if memq(car k,'(sin cos tan cot)) then
solvealgtrig(k,x) else
if memq(car k,'(sinh cosh tanh coth)) then
solvealghyp(k,x) else
if null x then t else
solvealggen(k,x)
) where x=solvealgtest(k);
symbolic procedure solvealgtest(k);
% Test if the arguments of a composite kernel interact with
% the variables known so far.
if atom k then nil else solvealgtest0(k);
symbolic procedure solvealgtest0(k);
% Test if kernel k interacts with the known variables.
solvealgtest1(k,iv!*) or solvealgtest1(k,uv!*);
symbolic procedure solvealgtest1(k,kl);
% list of those kernels in list kl, which occur somewhere
% in the composite kernel k.
if null kl then nil else
if member(k,kl) then list k else
if atom k then nil else
union(if smember(car kl,cdr k) then list car kl else nil,
solvealgtest1(k,cdr kl));
symbolic procedure solvealgvb k;
% Restricted variables are those which might establish
% non-algebraic relations like e.g. x + e**x. Test k
% and add it to the list.
fv!* := append(solvealgvb0 k,fv!*);
symbolic procedure solvealgvb0 k;
% test for restricted variables.
begin scalar ak;
ak := allkernels(k,nil);
if intersection(ak,iv!*) or intersection(ak,fv!*) then
error(99,list("transcendental variable dependency from",k));
return ak;
end;
symbolic procedure allkernels(a,kl);
% a is an algebraic expression. Extract all possible inner
% kernels of a and collect them in kl.
if numberp a then kl else
if atom a or a member uv!*
then if not member(a,kl) then a.kl else kl else
<<for each x in cdr a do
kl := allkernels1(numr s,allkernels1(denr s,kl))
where s=simp x;
kl>>;
symbolic procedure allkernels1(f,kl);
if domainp f then kl else
<<if not member(mvar f,kl) then
kl:=allkernels(mvar f,mvar f . kl);
allkernels1(lc f, allkernels1(red f,kl)) >>;
symbolic procedure solvealgexpt(k,x);
% kernel k is an exponential form.
( if eqcar(m,'quotient) and fixp caddr m then
if cadr m=1 then solvealgrad(cadr k,caddr m,x)
else solvealgradx(cadr k,cadr m,caddr m,x)
else if null x then solvealgid k
else if ((null intersection(w,uv!*) and
null intersection(w,iv!*) and
null intersection(w,fv!*))
where w=allkernels(m,nil))
then solvealggen(k,x)
else solvealgexptgen(k,x)
) where m = caddr k;
symbolic procedure solvealgexptgen(k,x);
% Kernel k is a general exponentiation u ** v.
begin scalar bas,xp,nv;
bas := cadr k; xp := caddr k;
if solvealgtest1(xp,uv!*) then return solvealgexptgen1(k,x)
else if solvealgtest1(bas,uv!*) then return solvealggen(k,x);
% remaining case: "constant" exponential expression to
% replaced by an id for syntatical reasons
nv := '(
% old kernel
( (expt !&alpha n))
% new variable
( !&beta)
% substitution
( ((expt !&alpha n) . !&beta) )
% inverse
( (!&beta (expt !&alpha n) !& ))
% new equations
nil
);
nv:=subst(bas,'!&alpha,nv);
nv:=subst(solve!-gensym(),'!&beta,nv);
nv:=subst(xp,'n,nv);
return solvealgupd(nv,nil);
end;
symbolic procedure solve!-gensym();
begin
scalar w;
w := explode solve!-gensymcounter;
solve!-gensymcounter := solve!-gensymcounter+1;
while length w < 4 do w := '!0 . w;
% If users have things to solve with names like G0001 in them, there
% could be confusion.
return compress ('g . w)
end;
symbolic procedure solvealgexptgen1(k,x);
% Kernel k is a general exponentiation u ** v.
% where v is an expression in a solution variable, u
% is constant. Transform all kernels with same basis
% and compatible exponent to common exponent denominator
% form.
begin scalar bas,xp,xpl,q,r,nk,sub;
bas := cadr k; xp := caddr k;
% collect all exponentials with this basis.
xpl:={(1 ./ 1).xp};
for each k in kl!* do
if eqcar(k,'expt) and cadr k=bas and
<<q:=simp{'quotient,r:=caddr k,xp};
fixp numr q and fixp denr q>> then
<<kl!*:=delete(k,kl!*); xpl:=(q.r).xpl>>;
% compute common denominator.
q:=1; for each e in xpl do q:=lcm(q,denr car e);
% the new artificial kernel.
nk:=reval{'expt,bas,{'quotient,xp,q}};
sub := for each e in xpl collect
{'expt,bas,cdr e}.
{'expt,nk,numr car e * q/denr car e};
system!*:=sublis(sub,system!*);
return solvealggen(nk,x);
end;
symbolic procedure solvealgradx(x,m,n,y);
% error(99,"forms e**(x/2) not yet implemented");
solvealgexptgen1({'expt,x,{'quotient,m,n}},y);
symbolic procedure solvealgrad(x,n,y);
% k is a radical exponentiation expression x**1/n.
begin scalar nv,m,!β
!&beta := solve!-gensym();
nv:= '(
% old kernel
( (expt !&alpha (quotient 1 !&n)))
% new variable
( !&beta)
% substitution
( ((expt !&alpha (quotient 1 !&n)) . !&beta) )
% inverse
% ( (!&beta !&alpha (expt !& !&n)) )
nil
% new equation
( (difference (expt !&beta !&n) !&alpha) )
);
m := list('!&alpha.x,'!&beta.!&beta,'!&n.n);
nv := subla(m,nv);
root!-vars!* := !&beta . root!-vars!*;
% prepare roots for simple surds.
if null y or y={x} then groebroots!* :=
({'plus,{'expt,!&beta,n},reval{'minus,x}}
.{{{'equal,!&beta,{'expt,x,{'quotient,1,n}}}}}).groebroots!*;
if null y then last!-vars!* := !&beta . last!-vars!*;
return solvealgupd(nv,y);
end;
symbolic procedure solvealgtrig0(f);
% examine if sin/cos identies must be applied.
begin scalar args,r,c;
args :=for each a in solvealgtrig01(f,nil) collect
(union(kernels numr q,kernels denr q) where q=simp a);
while args do
<<c:=car args;args:=cdr args;
for each q in args do r:=r or intersection(c,q)>>;
return r;
end;
symbolic procedure solvealgtrig01(f,args);
if atom f then args else
if memq(car f,'(sin cos tan cot sinh cosh tanh coth)) then
if constant_exprp cadr f then args else union({cadr f},args)
else solvealgtrig01(cdr f,solvealgtrig01(car f,args));
algebraic <<
operator p_sign,the_1;
let p_sign(~x) => if sign(x)=0 then 1 else sign(x);
let the_1(~x) =>1;
>>;
symbolic procedure solvealgtrig(k,x);
% k is a trigonometric function call.
begin scalar nv,m,s,!&alpha,!β
solvealgdb!* := union('(sin),solvealgdb!*);
if x then
if cdr x then
error(99,"too many variables in trig. function")
else x := car x;
solvealgvb k;
nv := '(
% old kernels
( (sin !&alpha)(cos !&alpha)(tan !&alpha)(cot !&alpha) )
% new variables
( (sin !&beta) (cos !&beta) )
% substitutions
( ((sin !&alpha) . (sin !&beta))
((cos !&alpha) . (cos !&beta))
%%% these should be handled now by the ruleset.
%%% ((tan !&alpha) . (quotient (sin !&beta) (cos !&beta)))
%%% ((cot !&alpha) . (quotient (cos !&beta) (sin !&beta)))
)
% inverses
(
((sin !&beta)
(cond ((and !*expli (test_trig))
'(!&loc (p_sign (!&!& !&))))
(t '(!&x (!&!& (root_of (equal (sin !&alpha)
!&) !&x))))))
((cos !&beta)
(cond ((and !*expli (test_trig))
'(!&x (plus (!&!& (times !&loc (acos !&)))
(times 2 pi !&arb))))
(t '(!&x (!&!& (root_of (equal (cos !&alpha)
!&) !&x))))))
)
% new equation
( (plus (expt (sin !&beta) 2)(expt (cos !&beta) 2) -1) )
);
% invert the inner expression.
s := if x then solvealginner(cadr k,x) else 'the_1;
!&beta := solve!-gensym();
m := list('!&alpha . (!&alpha:=cadr k),
'!&beta . !&beta,
'!&loc . solve!-gensym(),
'!&arb . {'arbint,!!arbint:=!!arbint+1},
'!&x . x,
'!&!& . s);
nv:=sublis!-pat(m , nv);
if x then last!-vars!*:=
append(last!-vars!*,{{'sin,!&beta},{'cos,!&beta}})
else const!-vars!* :=
append(const!-vars!*,{{'sin,!&beta}.{'sin,!&alpha},
{'cos,!&beta}.{'cos,!&alpha}});
return solvealgupd(nv,nil);
end;
symbolic procedure solvealghyp(k,x);
% k is a hyperbolic function call.
begin scalar nv,m,s,!&alpha,!β
solvealgdb!* := union('(sinh),solvealgdb!*);
if x then
if cdr x then
error(99,"too many variables in hyp. function")
else x := car x;
solvealgvb k;
nv := '(
% old kernels
( (sinh !&alpha)(cosh !&alpha)(tanh !&alpha)(coth !&alpha) )
% new variables
( (sinh !&beta) (cosh !&beta) )
% substitutions
( ((sinh !&alpha) . (sinh !&beta))
((cosh !&alpha) . (cosh !&beta))
)
% inverses
(
((sinh !&beta)
(cond ((and !*expli (test_hyp))
'(!&loc (p_sign (!&!& !&))))
(t '(!&x (!&!& (root_of (equal (sinh !&alpha)
!&) !&x))))))
((cosh !&beta)
(cond ((and !*expli (test_hyp))
'(!&x (plus (!&!& (times !&loc (acosh !&)))
(times 2 pi i !&arb))))
(t '(!&x (!&!& (root_of (equal (cosh !&alpha)
!&) !&x))))))
)
% new equation
( (plus (minus (expt (sinh !&beta) 2))(expt (cosh !&beta) 2) -1) )
);
% invert the inner expression.
s := if x then solvealginner(cadr k,x) else 'the_1;
!&beta := solve!-gensym();
m := list('!&alpha . (!&alpha:=cadr k),
'!&beta . !&beta,
'!&loc . solve!-gensym(),
'!&arb . {'arbint,!!arbint:=!!arbint+1},
'!&x . x,
'!&!& . s);
nv:=sublis!-pat(m , nv);
if x then last!-vars!*:=
append(last!-vars!*,{{'sinh,!&beta},{'cosh,!&beta}})
else const!-vars!* :=
append(const!-vars!*,{{'sinh,!&beta}.{'sinh,!&alpha},
{'cosh,!&beta}.{'cosh,!&alpha}});
return solvealgupd(nv,nil);
end;
symbolic procedure solvealgtrig2 u;
% r is a result from goesolve; remove trivial relations
% like sin^2 + cos^2 = 1.
begin scalar r,w,op,v,rh;
for each s in cdr u do
<<w := nil;
for each e in s do
% delete "sin u = sqrt(-cos u^2+1)" etc
if eqcar(e,'equal) and
(eqcar(cadr e,'sin) or eqcar(cadr e,'cos)) and
(op := caadr e) and (v := cadr cadr e) and
member(if eqcar(rh:=caddr e,'!*sq!*) then cadr rh else rh,
subst({if op='sin then 'cos else 'sin,v},'!-form!-,
'((MINUS (SQRT (PLUS (MINUS (EXPT !-form!- 2)) 1)))
(SQRT (PLUS (MINUS (EXPT !-form!- 2)) 1)))))
then nil
else w:=e.w;
w := reverse w;
if not member(w,r) then r:=w.r;
>>;
return 'list . reverse r;
end;
symbolic procedure solvealghyp2 u;
% r is a result from goesolve; remove trivial relations
% like cosh^2 - sinh^2 = 1.
begin scalar r,w,op,v,rh;
for each s in cdr u do
<<w := nil;
for each e in s do
% delete "sinh u = sqrt(cosh u^2-1)","cosh u = sqrt(sinh u^2+1)"
if eqcar(e,'equal) and
(eqcar(cadr e,'sinh) or eqcar(cadr e,'cosh)) and
(op := caadr e) and (v := cadr cadr e) and
member(if eqcar(rh:=caddr e,'!*sq!*) then cadr rh else rh,
if op='sinh then
subst({'cosh,v},'!-form!-,
'((MINUS (SQRT (PLUS (EXPT !-form!- 2) 1)))
(SQRT (PLUS (EXPT !-form!- 2) 1))))
else
subst({'sinh,v},'!-form!-,
'((MINUS (SQRT (PLUS (EXPT !-form!- 2) (MINUS 1))))
(SQRT (PLUS (EXPT !-form!- 2) (MINUS 1))))))
then nil
else w:=e.w;
w := reverse w;
if not member(w,r) then r:=w.r;
>>;
return 'list . reverse r;
end;
symbolic procedure solvealggen(k,x);
% k is a general function call; processable if SOLVE
% can invert the function.
begin scalar nv,m,s;
if cdr x then
error(99,"too many variables in function expression");
x := car x;
solvealgvb k;
nv := '(
% old kernels
( !&alpha )
% new variables
( !&beta )
% substitutions
( ( !&alpha . !&beta) )
% inverses
(( !&beta '(!&x (!&!& !&))))
% new equation
nil);
% invert the kernel expression.
s := solvealginner(k,x);
m := list('!&alpha . k,
'!&beta . solve!-gensym(),
'!&x . x,
'!&!& . s);
nv:=sublis!-pat(m , nv);
return solvealgupd(nv,nil);
end;
symbolic procedure solvealgid k;
% k is a "constant" kernel, however in a syntax unprocessable
% for Groebner (e.g. expt(a/2)); replace temporarily
begin scalar nv,m;
nv := '(
% old kernels
( !&alpha )
% new variables
( )
% substitutions
( ( !&alpha . !&beta) )
% inverses
(( !&beta nil . !&alpha))
% new equation
nil);
% invert the kernel expression.
m := list('!&alpha . k, '!&beta . solve!-gensym());
nv:=sublis(m , nv);
return solvealgupd(nv,nil);
end;
symbolic procedure solvealginner(s,x);
<<s := solveeval1 {{'equal,s,'!#},{'LIST,x}};
s := reval cadr s;
if not eqcar(s,'EQUAL) or not equal(cadr s,x) then
error (99,"inner expression cannot be inverted");
{'lambda,'(!#),caddr s}>>;
symbolic procedure solvealgupd(u,innervars);
% Update the system and the structures.
begin scalar ov,nv,sub,inv,neqs;
ov := car u; u := cdr u;
nv := car u; u := cdr u;
sub:= car u; u := cdr u;
inv:= car u; u := cdr u;
neqs:=car u; u := cdr u;
for each x in ov do kl!*:=delete(x,kl!*);
for each x in innervars do
for each y in nv do depend1(y,x,t);
sub!* := append(sub,sub!*);
iv!* := append(nv,iv!*);
inv!* := append(inv,inv!*);
system!* := append(
for each u in neqs collect
<<u:= numr simp u; solvealgk0 u; u>>,
for each u in system!* collect numr subf(u,sub) );
return t;
end;
symbolic procedure solvealginv u;
% Reestablish the original variables, produce inverse
% mapping and do complete value propagation.
begin scalar v,r,s,m,lh,rh,y,z,tag,sub0,sub,!*expli,noarb,arbs;
scalar abort;
integer n;
sub0 := for each p in sub!* collect (cdr p.car p);
tag := t;
r := for each sol in cdr u join
<<sub := sub0; abort := v:= r:= s:= noarb :=arbs :=nil;
if !*test_solvealg then
<<prin2t "================================";
prin2t const!-vars!*;
prin2t " next basis:";
writepri(mkquote sol,'only);
>>;
for each eqn in reverse cdr sol do
<<lh := cadr eqn; rh := subsq(simp!* caddr eqn,s);
if !*test_solvealg then
writepri(mkquote {'equal,lh,prepsq rh},'only);
!*expli:=member(lh,iv!*);
% look for violated constant relations.
if (y:=assoc(lh,const!-vars!*)) and constant_exprp prepsq rh
and numr subtrsq(rh,simp cdr y) then abort := t;
% look for a "negative" root.
if memq(lh,root!-vars!*) and numberp(y:=reval{'sign,prepsq rh})
and y<0 then abort := t;
if not !*expli then noarb := t;
if !*expli and not noarb then
<< % assign value to free variables;
for each x in uv!* do
if !*arbvars and solvealgdepends(rh,x) and not member(x,fv!*)
and not member(x,arbs) then
<<z := mvar makearbcomplex();
y := z; v := x . v; r := simp y . r;
% rh := subsq(rh,list(x.y));
% s := (x . y) . s;
arbs:=x.arbs;
>>;
if not smemq('root_of,rh) then
s:=(lh.prepsq rh).s
else fv!*:=lh.fv!*;
>>;
if (m:=assoc(lh,inv!*))then
<<m:=cdr m; lh :=car m; kl!* := eqn;
if eqcar(lh,'cond) or eqcar(lh,'quote) then
lh:=car(m:=eval lh);
rh:=solvenlnrsimp subst(prepsq rh,'!&,cadr m)>>;
% if local variable, append to substitution.
if not member(lh,uv!*) and !*expli then
<< sub:=append(sub,{lh .(z:=prepsq subsq(rh,sub))});
if smember(lh,r) then r:=subst(z,lh,r)
>>;
% append to the final output.
if (member(lh,uv!*) or not !*expli)
% inhibit repeated same values.
and not<< z:=subsq(rh,sub);
n:=length member(z,r);
n>0 and lh=nth(v,length v + 1 - n)>>
then <<r:=z.r; v:=lh.v;>>;
>>;
% Classify result.
% for each x in uv!* do
% if tag and not member(x,v) and smember(x,r) then tag:=nil;
if !*test_solvealg then
if abort then yesp "ABORTED" else
<<prin2t " --------> ";
writepri(mkquote ('list .for each u in pair(v,r) collect
{'equal,car u,prepsq cdr u}) ,'only);
prin2t "================================";
yesp "continue?";
>>;
if not abort then {reverse r . reverse v}
>>;
return solvealg!-verify(tag,r);
end;
symbolic procedure solvealgdepends(u,x);
% inspect u for explicit dependency of x, being careful for
% root_of subexpressions.
if u=x then t else if atom u then nil else
if eqcar(u,'root_of) then
if x=caddr u then nil else solvealgdepends(cadr u,x) else
solvealgdepends(car u,x) or solvealgdepends(cdr u,x);
symbolic procedure test_trig();
begin scalar lh,rh,r;
lh := cadr kl!*; rh:= caddr kl!*;
if member(lh . nil, solvealgdb!*) then return nil;
r := not !*complex and not smemq('i,kl!*) and
not smemq('!:gi!:,kl!*) and not smemq('!:cr!:,kl!*) and
not smemq('root_of,kl!*);
if not r then solvealgdb!* :=
append(solvealgdb!*,{('sin.cdr lh).nil,('cos.cdr lh).nil});
return r;
end;
symbolic procedure test_hyp();
begin scalar lh,rh,r;
lh := cadr kl!*; rh:= caddr kl!*;
if member(lh . nil, solvealgdb!*) then return nil;
r := not !*complex and not smemq('i,kl!*) and
not smemq('!:gi!:,kl!*) and not smemq('!:cr!:,kl!*) and
not smemq('root_of,kl!*);
if not r then solvealgdb!* :=
append(solvealgdb!*,{('sinh.cdr lh).nil,('cosh.cdr lh).nil});
return r;
end;
fluid '(!*solvealg_verify);
% the idea of the following procedure is to exclude isolated
% solutions which give a substantial residue when subsituted
% into the equation system under "on rounded"; as long as no
% good criterion for a residue to be small has been found, this
% step is disabled.
symbolic procedure solvealg!-verify(tag,r);
<<if !*rounded and !*solvealg_verify then
begin scalar min,s,cmpl,!*msg;
% exclude solutions with a residue substantially
% above the minimum of all nonzero residues.
cmpl:=!*complex;
if not cmpl then setdmode('complex,!*complex:=t);
s := for each u in r collect solvealg!-verify1 u.u;
min:=simp'(quotient 1 100);
r:= for each u in s join
if null car u or minusf numr subtrsq(car u,min) then {cdr u};
if not cmpl then
<<setdmode('complex,nil); !*complex:=nil>>;
end;
tag . for each q in r collect car q . cdr q . list 1
>>;
symbolic procedure solvealg!-verify1 s;
% verify solution s for the current equation system.
begin scalar sub,nexpli,x,y,sum,fail;
sub:= for each u in pair(cdr s,car s) collect if not nexpli then
<<y:=prepsq cdr u;
if not (domainp y or constant_exprp y) then nexpli:=t;
car u.y>>;
% a non explicit solution cannot be tested.
if nexpli then return nil;
sum := nil ./ 1;
for each u in osystem!* do if not fail then
<<x:=subf(u,sub);
if domainp numr x then
sum:=addsq(sum,absf numr x ./ denr x)
else fail := t;
>>;
return if fail then nil else sum;
end;
symbolic procedure sublis!-pat(a,u);
% like sublis, but replace lambda expressions by matching their
% actual arguments.
begin scalar v;
if atom u then return
<<v:=assoc(u,a); if v then sublis!-pat(a,cdr v) else u>>;
v:=assoc(car u,a);
if v and (v:=cdr v) and eqcar(v,'lambda) then return
sublis!-pat((caadr v.cadr u).a,caddr v);
return sublis!-pat1(a,u);
end;
symbolic procedure sublis!-pat1(a,l);
if null l then nil else
if atom l then sublis!-pat(a,l) else
sublis!-pat(a,car l) . sublis!-pat1(a,cdr l);
%----------------------------------------------------------------
% section for single trigonometric polynomials
%----------------------------------------------------------------
symbolic procedure solvenonlnrtansub(p,x);
% Perform tangent substitution.
if not smemq('sin,p) and not smemq('cos,p) then
if smemq(x,p) then nil else nil.p
else if car p='cos then
if smemq(x,cdr p) then (cdr p).
'(quotient (difference 1(expt tg!- 2)) (plus 1(expt tg!- 2)))
else nil.p
else if car p='sin then
if smemq(x,cdr p) then (cdr p).
'(quotient (times 2 tg!-) (plus 1(expt tg!- 2)))
else nil.p
else
(if ca and cd and
(car ca = car cd or null car ca or null car cd)
then (car ca or car cd).(cdr ca.cdr cd))
where ca=solvenonlnrtansub(car p,x),
cd=solvenonlnrtansub(cdr p,x);
symbolic procedure solvenonlnrtansolve(u,x,w);
begin scalar v,s,z,r,y;
integer ar;
% We reset arbint for each solve call such that equal forms can
% be recognized by the function union.
ar := !!arbint;
v:=caar u; u:= prepf numr simp cdr u;
s:=solveeval {u,'tg!-};
!!arbint:=ar;
for each q in cdr s do
<<z:=reval caddr q;
z:=reval sublis(solvenonlnrtansolve1 z,z);
!!arbint:=ar;
y:=solve0({'equal,{'tan,{'quotient,V,2}},z},x);
r:=union(y,r);
>>;
% test for the special cases x=pi(not covered
% by tangent substitution).
if null numr subf(w,{x.'pi})
then <<!!arbint:=ar; r:=union(solve0({'equal,{'cos,x},-1},x),r)>>;
return t.r;
end;
symbolic procedure solvenonlnrtansolve1 u;
% Find all cos**2.
if atom u then nil else
if car u='expt and eqcar(cadr u,'cos) and caddr u=2 then
{u . {'difference,1,{'expt,{'sin,cadr cadr u},2}}}
else union(solvenonlnrtansolve1 car u,solvenonlnrtansolve1 cdr u);
%----------------------------------------------------------------
% section for single hyperbolic polynomials
%----------------------------------------------------------------
symbolic procedure solvenonlnrtanhsub(p,x);
% Perform hyperbolic tangent substitution.
if not smemq('sinh,p) and not smemq('cosh,p) then
if smemq(x,p) then nil else nil.p
else if car p='cosh then
if smemq(x,cdr p) then (cdr p).
'(quotient (plus 1(expt tgh!- 2)) (difference 1(expt tgh!- 2)))
else nil.p
else if car p='sinh then
if smemq(x,cdr p) then (cdr p).
'(quotient (times 2 tgh!-) (difference 1(expt tgh!- 2)))
else nil.p
else
(if ca and cd and
(car ca = car cd or null car ca or null car cd)
then (car ca or car cd).(cdr ca.cdr cd))
where ca=solvenonlnrtanhsub(car p,x),
cd=solvenonlnrtanhsub(cdr p,x);
symbolic procedure solvenonlnrtanhsolve(u,x,w);
begin scalar v,s,z,r,y,ar;
ar := !!arbint;
v:=caar u; u:= prepf numr simp cdr u;
s:=solveeval {u,'tgh!-};
ar := !!arbint;
for each q in cdr s do
<<z:=reval caddr q;
z:=reval sublis(solvenonlnrtanhsolve1 z,z);
!!arbint:=ar;
y:=solve0({'equal,{'tanh,{'quotient,v,2}},z},x);
r:=union(y,r);
>>;
if !*complex and null numr subf(w,{x.'(times pi i)})
then <<!!arbint:=ar; r:=union(solve0({'equal,{'cosh,x},-1},x),r)>>;
return t.r;
end;
symbolic procedure solvenonlnrtanhsolve1 u;
% Find all cosh**2.
if atom u then nil else
if car u='expt and eqcar(cadr u,'cosh) and caddr u=2 then
{u . {'plus,1,{'expt,{'sinh,cadr cadr u},2}}}
else union(solvenonlnrtanhsolve1 car u,solvenonlnrtanhsolve1 cdr u);
endmodule;
end;