Artifact 2bb70dd3eda88bcc518e4aea94f7ac3ba3688971177c7af6d1a79a171d586a44:
- Executable file
r36/src/solve.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 111400) [annotate] [blame] [check-ins using] [more...]
module /ure solsolve; % Solve one or more algebraic equations. % Author: David R. Stoutemyer. % Major modifications by: David Hartley, Anthony C. Hearn, Herbert % Melenk, Donald R. Morrison and Rainer Schoepf. create!-package('(solve ppsoln solvelnr glsolve solvealg solvetab quartic),nil); % Other packages needed by solve package. load!-package 'matrix; fluid '(!*allbranch !*arbvars !*exp !*ezgcd !*fullroots !*limitedfactors !*multiplicities !*notseparate !*numval !*numval!* !*rounded !*solvealgp !*solvesingular !*varopt !!gcd !:prec!: asymplis!* alglist!* dmode!* kord!* vars!* !*!*norootvarrenamep!*!*); % NB: !*!*norootvarrenamep!*!* is internal to this module, and should % *never* be changed by a user. fluid '(inside!-solveeval solve!-gensymcounter); solve!-gensymcounter := 1; global '(!!arbint multiplicities!* assumptions requirements); switch allbranch,arbvars,fullroots,multiplicities,solvesingular; % varopt,nonlnr. !*varopt := t; put('fullroots,'simpfg,'((t (rmsubs)))); flag('(!*allbranch multiplicities!* assumptions requirements), 'share); % ***** Some Non-local variables ***** !*allbranch := t; % Returns all branches of solutions if T. !*arbvars := t; % Presents solutions to singular systems % in terms of original variables if NIL % !*multiplicities Lists all roots with multiplicities if on. % !*fullroots := t; % Computes full roots of cubics and quartics. !*solvesingular := t; % Default value. % !!gcd SOLVECOEFF returns GCD of powers of its arg in % this. With the decompose code, this should % only occur with expressions of form x^n + c. algebraic operator one_of; put('arbint,'simpfn,'simpiden); % algebraic operator arbreal; symbolic operator expand_cases; symbolic procedure simp!-arbcomplex u; simpiden('arbcomplex . u) where dmode!*=nil; deflist('((arbcomplex simp!-arbcomplex)),'simpfn); % ***** Utility Functions ***** symbolic procedure freeofl(u,v); null v or freeof(u,car v) and freeofl(u,cdr v); symbolic procedure allkern elst; % Returns list of all top-level kernels in the list of standard % quotients elst. Corrected 5 Feb 92 by Francis Wright. if null elst then nil else union(kernels numr car 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. We need to reorder coefficients % since kernel order can change in the calling function. 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 . reorder lc ex) . ans; ex := red ex>>; if ex then ans := (0 . reorder ex) . ans; return ans end; % ***** Evaluation Interface ***** % Solvemethods!* is a list of procedures which are able to process % one problem class. Each of its members must check itself % whether it can be applied or not. The classical equation solver % is called if none of the methods can contribute. % % Protocol: % % input: PSOPFN standard, where the elements of the input list % have been passed through REVAL. % % output: % 'nil: the algorithm cannot be applied because the problem % belongs to a different problem class; % '(failed): the problem belongs to the class represented % by the algorithm but the program has been % unable to compute a result. The problem should % not be given to any other method - instead the % input should be returned. % result: the algorithm has been successful and the final % result is returned as algebraic form (including an % eventually empty result for an "inconsistent" case). fluid '(solvemethods!*); put('solve,'psopfn,'solvemaster); symbolic procedure solvemaster u; begin scalar w,r,m; w:=for each q in u collect reval q; m:=solvemethods!*; while null r and m do <<r:=apply1(car m,w); m:=cdr m>>; return if null r then solveeval w else if eqcar(r,'failed) then 'solve.u else r; end; symbolic procedure solveeval u; begin scalar !*ezgcd,!!gcd,vars!*; integer nargs; if atom u then rerror(solve,1,"SOLVE called with no equations") else if null dmode!* then !*ezgcd := t; nargs := length u; if not inside!-solveeval then <<solve!-gensymcounter := 1; assumptions :=requirements:={'list}>>; u := (if nargs=1 then solve0(car u,nil) else if nargs=2 then solve0(car u, cadr u) else <<lprim "Please put SOLVE unknowns in a list"; solve0(car u,'list . cdr u)>>) where inside!-solveeval = t, !*resimp = not !*exp; if not inside!-solveeval then <<assumptions := solve!-clean!-info(assumptions,t); requirements:= solve!-clean!-info(requirements,nil)>>; return !*solvelist2solveeqlist u 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 !*solvelist2solveeqlist u; begin scalar x,y,z; for each j in u do <<if caddr j=0 then rerror(solve,2,"zero multiplicity") else if null cadr j then x := for each k in car j collect list('equal,!*q2a k,0) else x := for each k in pair(cadr j,car j) collect list('equal,car k,!*q2a cdr k); if length vars!* > 1 then x := 'list . x else x := car x; z := (caddr j . x) . z>>; z := sort(z,function ordp); x := nil; if !*multiplicities then <<for each k in z do for i := 1:car k do x := cdr k . x; multiplicities!* := nil>> else <<for each k in z do << x := cdr k . x; y := car k . y>>; multiplicities!* := 'list . reversip y>>; % Now check for redundant solutions. % if length vars!*>1 then z := check_solve_redundancy z; return 'list . reversip x end; % symbolic procedure check_solve_redundancy u; % % We assume all solutions are prefixed by LIST. % begin scalar x,y; % x := for each j in u collect cdr j; % Remove the LIST. % for each j in u do if not supersetlist(cdr j,x) then y:= j . y; % return reversip!* y % end; symbolic procedure supersetlist(u,v); % Returns true if u is a non-equal superset of any element of v. v and (u neq car v and null setdiff(car v,u) or supersetlist(u,cdr v)); % ***** Fundamental SOLVE Procedures ***** Comment most of these procedures return a list of "solve solutions". A solve solution is a list with three fields: the list of solutions, the corresponding variables (or NIL if the equations could not be solved --- in which case there is only one solution in the first field) and the multiplicity; symbolic procedure solve0(elst,xlst); % This is the driving function for the solve package. % Elst is any prefix expression, including a list prefixed by LIST. % Xlst is a kernel or list of kernels. Solves eqns in elst for % vars in xlst, returning either a list of solutions, a single % solution, or NIL if the solutions are inconsistent. begin scalar !*exp,!*notseparate,w; integer neqn; !*exp := !*notseparate := t; % Form a list of equations as expressions. elst := for each j in solveargchk elst collect simp!* !*eqn2a j; neqn := length elst; % There must be at least one. % Determine variables. 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; % More than one equation or variable, or single eqn has no solution. 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; symbolic procedure basic!-kern u; <<for each k in u do w:=union(basic!-kern1 k,w); w>> where w=nil; symbolic procedure basic!-kern1 u; % Expand a composite kernel. begin scalar w; if atom u then return {u} else if algebraic!-function car u and (w := allbkern for each q in cdr u collect simp q) then return w else return {u} end; symbolic procedure algebraic!-function q; % Returns T if q is a name of an operator with algebraic evaluation % properties. flagp(q,'realvalued) or flagp(q,'alwaysrealvalued) or get(q,'!:rd!:) or get(q,'!:cr!:) or get(q,'opmtch); symbolic procedure allbkern elst; % extract all elementary kernels from list of quotients. if null elst then nil else union(basic!-kern kernels numr car elst, allbkern cdr elst); symbolic procedure solvevars elst; <<for each j in allbkern elst do if not constant_exprp j then s := ordad(j,s); s>> where s=nil; symbolic procedure solutionp u; null u or cadar u and not root_of_soln_p caar u and solutionp cdr u; symbolic procedure root_of_soln_p u; null cdr u and kernp (u := car u) and eqcar(mvar numr u,'root_of); symbolic procedure solveargchk u; if getrtype (u := reval 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 collecting side information about the solution. symbolic procedure solve!-clean!-info(fl,flg); % Check for constants and multiples in side relations fl. % If flg is t then relations are factorised and constants removed. % Otherwise the relations are autoreduced and the presence of a % constant truncates the whole list. 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:= sqfrf!* absf numr simp{'nprimitive,form}; 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; % ***** Procedures for solving a single eqn ***** %symbolic procedure solvesq (ex,var,mul); % begin scalar r; % r:= for each w in solvesq1(ex,var,mul) join % if null cadr w % or numr subf(denr ex,{caadr w.prepsq caar w}) then {w}; % if r and not domainp denr ex then % assumptions:=append(assumptions,{prepf denr ex}); % return r; % end; % Patch from 12 Feb 05 symbolic procedure solvesq (ex,var,mul); begin scalar r,x; r:= for each w in solvesq1(ex,var,mul) join if null cadr w or eqcar(x := prepsq caar w,'root_of) or numr subfx(denr ex,{caadr w . x}) then {w}; if r and not domainp denr ex then assumptions:=append(assumptions,{prepf denr ex}); return r end; symbolic procedure solvesq1 (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,oldkorder,x1,y,z; integer mu; ex := numr ex; if null(x1 := topkern(ex,var)) then return nil; oldkorder := setkorder list var; % The following section should be extended for other cases. 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}}; % don't call fctrf here in rounded mode, so polynomial won't get % rounded (since factoring isn't going to succeed) ex := if !*rounded then {1,ex . 1} else 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; % Test for decomposition of e1. Only do if rounded is off. 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 smemq(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 := append( 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 list list(list(e1 ./ 1),nil,mu), else mkrootsof(e1 ./ 1,var,mu), z)>>; return z end; symbolic procedure solvedecomp(u,var,mu); % Solve for decomposed expression. At the moment, only one % level of decomposition is considered. begin scalar failed,x; if length(x := solve0(car u,cadadr u))=1 then return nil; u := cdr u; while u do <<x := for each j in x conc if caddr j neq 1 or null cadr j then <<lprim list("Tell Hearn solvedecomp",x,u); failed := t; nil>> else solve0(list('difference,prepsq caar j,caddar u), if cdr u then cadadr u else var); if failed then u := nil else u := cdr u>>; return if failed then nil else adjustmul(x,mu) end; symbolic procedure adjustmul(u,n); % Multiply the multiplicities of the solutions in u by n. if n=1 then u else for each x in u collect list(car x,cadr x,n*caddr x); 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 a list of solve solutions; begin scalar !*numval!*; !*numval!* := !*numval; % Keep value for use in solve11. return solve11(e1,x1,var,mu) end; symbolic procedure solve11(e1,x1,var,mu); begin scalar !*numval,b,coefs,hipow,w; integer n; % The next test should check for true degree in var. if null !*fullroots and null !*rounded and numrdeg(e1,var)>2 then return mkrootsof(e1 ./ 1,var,mu); !*numval := t; % Assume that actual numerical values wanted. coefs:= errorset!*(list('solvecoeff,mkquote e1,mkquote x1),nil); % if atom coefs then return list list(list(e1 . 1),nil,mu); if atom coefs or atom x1 and x1 neq var then return mkrootsof(e1 ./ 1,var,mu); % solvecoeff problem - no soln. coefs := car coefs; n:= !!gcd; % numerical gcd of powers. 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)); % 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 := 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)) % 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 solvequadratic(getcoeff(coefs,2), getcoeff(coefs,1),getcoeff(coefs,0)) conc solvesq(subtrsq(x1,j),var,mu)>> else return solvehipow(e1,x1,var,mu,coefs,hipow) end; symbolic procedure getcoeff(u,n); % Get the nth coefficient in the list u as a standard quotient. if null u then nil ./ 1 else if n=caar u then cdar u else if n<caar u then nil ./ 1 else getcoeff(cdr u,n); symbolic procedure putcoeff(u,n,v); % Replace the nth coefficient in the list u by v. if null u then list(n . v) else if n=caar u then (n . v) . cdr u else if n<caar u then (n . v) . u else car u . putcoeff(cdr u,n,v); symbolic procedure solvehipow(e1,x1,var,mu,coefs,hipow); % Solve a system with degree greater than 2. Since we cannot write % down the solution directly, we look for various forms that we % know how to solve. begin scalar b,c,d,f,rcoeffs; f:=(hipow+1)/2; d:=exptsq(simp!* x1,!!gcd); rcoeffs := reverse coefs; return if solve1test1(coefs,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 <<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) conc solvesq(addsq(1 ./ 1,multsq(d,subtrsq(d,caar j))), var,caddr j)>> else if solve1test2(coefs,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 % now look for real roots before cubics or quartics. We must % reverse the answer from solveroots so that roots come out % in same order from SOLVE. % else if !*numval!* and dmode!* memq '(!:rd!: !:cr!:) % this forces solveroots independent of numval. 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)) conc solvesq(subtrsq(d,j),var,mu) 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)) conc solvesq(subtrsq(d,j),var,mu) % else list list(list(e1 ./ 1),nil,mu) else mkrootsof(e1 ./ 1,var,mu) % We can't solve quintic and higher. 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(coefs,rcoeffs,f); % True if equation is symmetric in its coefficients. f is midpoint. begin integer j,p; if null coefs or caar coefs neq 0 or null !*fullroots then return nil; p := caar coefs + caar rcoeffs; a: if j>f then return t else if (caar coefs + caar rcoeffs) neq p or cdar coefs neq cdar rcoeffs then return nil; coefs := cdr coefs; rcoeffs := cdr rcoeffs; j := j+1; go to a end; symbolic procedure solve1test2(coefs,rcoeffs,f); % True if equation is antisymmetric in its coefficients. f is % midpoint. begin integer j,p; if null coefs or caar coefs neq 0 or null !*fullroots then return nil; p := caar coefs + caar rcoeffs; a: if j>f then return t else if (caar coefs + caar rcoeffs) neq p or numr addsq(cdar coefs,cdar rcoeffs) then return nil; coefs := cdr coefs; 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; % the following line made solve(x^(1/2)=0) etc. wrong % if null numr lincoeff then return nil; u := car u; return if freeof(car u,var) % c**(...) = b. then if null numr lincoeff then nil else <<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(solveexpt!-logterm lincoeff, simp!* c), simp!* list('log,car u))),var,mu)>> else if freeof(cadr u,var) and null numr lincoeff %(...)**b=0. then if check!-condition {'equal,{'sign,cadr u},1} then solvesq(simp!* car u,var,mu) else solveexpt!-rootsof(u,lincoeff,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('quotient, list('times,2,'pi, list('arbint,!!arbint)), cadr u)>> % 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) else solveexpt!-rootsof(u,lincoeff,var,mu) end; symbolic procedure solveexpt!-rootsof(u,lincoeff,var,mu); mkrootsof(subtrsq(simp!*('expt . u),lincoeff),var,mu); put('expt,'solvefn,'solveexpt); symbolic procedure solveexpt!-logterm lincoeff; % compute log(lincoeff), ignoring multiplicity and converting % log(-a) to log(a) + i pi. simp!* list('log,mk!*sq lincoeff); % if not !*allbranch or not minusf numr lincoeff % then simp!* list('log,mk!*sq lincoeff) % else % addsq(simp!*'(times i pi), % simp!* {'log,mk!*sq(negf numr lincoeff ./ denr lincoeff)}); 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 solveinvpat(u,op); begin scalar c,f; f:=get(op,'solveinvpat); if smemq('arbint,f) then f:=subst( if !*allbranch then list('arbint,!!arbint:=!!arbint+1) else 0, 'arbint,f); if not !*allbranch then f:={car f}; return for each c in reverse f join solvesq(simp!* subst(caar u,'(~v),subst(mk!*sq cadddr u,'(~r),c)), cadr u,caddr u) end; put('cos,'solveinvpat, { quote (- ~v + acos(~r) + 2*arbint*pi), quote (- ~v - acos(~r) + 2*arbint*pi) }); put('cos,'solvefn, '(lambda(u) (solveinvpat u 'cos))); put('sin,'solveinvpat, { quote (- ~v + asin(~r) + 2*arbint*pi), quote (- ~v - asin(~r) + 2*arbint*pi + pi) }); put('sin,'solvefn, '(lambda(u) (solveinvpat u 'sin))); put('sec,'solveinvpat, { quote (- ~v + asec(~r) + 2*arbint*pi), quote (- ~v - asec(~r) + 2*arbint*pi) }); put('sec,'solvefn, '(lambda(u) (solveinvpat u 'sec))); put('csc,'solveinvpat, { quote (- ~v + acsc(~r) + 2*arbint*pi), quote (- ~v - acsc(~r) + 2*arbint*pi + pi) }); put('csc,'solvefn, '(lambda(u) (solveinvpat u 'csc))); put('tan,'solveinvpat, { quote (- ~v + atan(~r) + arbint*pi)}); put('tan,'solvefn, '(lambda(u) (solveinvpat u 'tan))); put('cot,'solveinvpat, { quote (- ~v + acot(~r) + arbint*pi)}); put('cot,'solvefn, '(lambda(u) (solveinvpat u 'cot))); put('cosh,'solveinvpat, { quote (- ~v + acosh(~r) + 2*arbint*i*pi), quote (- ~v - acosh(~r) + 2*arbint*i*pi) }); put('cosh,'solvefn, '(lambda(u) (solveinvpat u 'cosh))); put('sinh,'solveinvpat, { quote (- ~v + asinh(~r) + 2*arbint*i*pi), quote (- ~v - asinh(~r) + 2*arbint*i*pi + i*pi) }); put('sinh,'solvefn, '(lambda(u) (solveinvpat u 'sinh))); put('sech,'solveinvpat, { quote (- ~v + asech(~r) + 2*arbint*i*pi), quote (- ~v - asech(~r) + 2*arbint*i*pi) }); put('sech,'solvefn, '(lambda(u) (solveinvpat u 'sech))); put('csch,'solveinvpat, { quote (- ~v + acsch(~r) + 2*arbint*i*pi), quote (- ~v - acsch(~r) + 2*arbint*i*pi + i*pi) }); put('csch,'solvefn, '(lambda(u) (solveinvpat u 'csch))); put('tanh,'solveinvpat, { quote (- ~v + atanh(~r) + arbint*i*pi)}); put('tanh,'solvefn, '(lambda(u) (solveinvpat u 'tanh))); put('coth,'solveinvpat, { quote (- ~v + acoth(~r) + arbint*i*pi)}); put('coth,'solvefn, '(lambda(u) (solveinvpat u 'coth))); symbolic procedure mkexp u; reval(aeval!*({'plus,{'cos,x},{'times,'i,{'sin,x}}} where x = reval u) where !*rounded = nil,dmode!* = nil); symbolic procedure solvecoeff(ex,var); % Ex is a standard form and var a kernel. Returns a list of % dotted pairs of exponents and coefficients (as standard quotients) % of var in ex, lowest power first, with exponents divided by their % gcd. This gcd is stored in !!GCD. begin scalar clist,oldkord; oldkord := updkorder var; clist := reorder ex; setkorder oldkord; clist := coeflis clist; !!gcd := caar clist; for each x in cdr clist do !!gcd := gcdn(car x,!!gcd); for each x in clist do <<rplaca(x,car x/!!gcd); rplacd(x,cdr x ./ 1)>>; return clist end; symbolic procedure solveroots(ex,var,mu); % Ex is a square and content free univariate standard form, var the % relevant variable and mu the root multiplicity. Finds insoluble, % complex roots of EX, returning a list of solve solutions. begin scalar y; y := reval list('root_val,mk!*sq(ex ./ 1)); if not(car y eq 'list) then errach list("incorrect root format",ex); return for each z in cdr y collect % if not(car z eq 'equal) or cadr z neq var % then errach list("incorrect root format",ex) % else list(list simp caddr z,list var,mu) list(list simp z,list var,mu) end; % ***** Procedures for solving a system of eqns ***** Comment. The routines for solving systems of equations return a "tagged solution list", where tagged solution list ::= tag . list of solve solution tag ::= t | nil | 'inconsistent | 'singular | 'failed solve solution ::= {solution rhs,solution lhs,multiplicity} | {solution rhs,nil,multiplicity} solution rhs ::= list of sq solution lhs ::= list of kernel multiplicity ::= posint If the tag is anything but t, the list of solve solutions is empty. See solvenonlnrsys for more about the tags; symbolic procedure solvesys(exlis,varlis); % exlis: list of sf, varlis: list of kernel % -> solvesys: tagged solution list % The expressions in exlis are reordered wrt the kernels in varlis, % and solved. For some switch settings, the internal % solve procedure may produce an error, leaving the kernel order % disturbed, so an errorset is used here. begin scalar oldkord; % The standard methods for linear and polynomial system % don't work for non-prime modulus. if !*modular then <<load!-package 'modsr; return msolvesys(exlis,varlis,t)>>; oldkord := setkorder varlis; exlis := for each j in exlis collect reorder j; exlis := errorset!*({'solvemixedsys,mkquote exlis,mkquote varlis}, t); setkorder oldkord; if errorp exlis then error1(); return car exlis; end; symbolic procedure solvemixedsys(exlis,varlis); % exlis: list of sf, varlis: list of kernel % -> solvemixedsys: tagged solution list % Solve a mixed linear/nonlinear system, solving the linear % part and substituting into the nonlinear until a core nonlinear % system remains. Assumes solvenonlnrsys and solvelnrsys both handle % all trivial cases properly. if null cadr(exlis := siftnonlnr(exlis,varlis)) then % linear solvelnrsys(car exlis,varlis) else if null car exlis then % nonlinear solvenonlnrsys(cadr exlis,varlis) else % mixed 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 subf(ex,z) then {ex}; % resimp?? varlis := setdiff(varlis,cadr x); % remaining free variables 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; symbolic procedure siftnonlnr(exlis,varlis); % exlis: list of sf, varlis: list of kernel % -> siftnonlnr: {list of sf, list of sf} % separate exlis into {linear,nonlinear} begin scalar lin,nonlin; foreach ex in exlis do if nonlnr(ex,varlis) then nonlin := ex . nonlin else lin := ex . lin; return {reversip lin,reversip nonlin}; end; symbolic procedure nonlnrsys(exlis,varlis); % exlis: list of sf, varlis: list of kernel % -> nonlnrsys: bool if null exlis then nil else nonlnr(car exlis,varlis) or nonlnrsys(cdr exlis,varlis); symbolic procedure nonlnr(ex,varlis); % ex: sf, varlis: list of kernel -> nonlnr: bool if domainp ex then nil else if mvar ex member varlis then ldeg ex>1 or not freeofl(lc ex,varlis) or nonlnr(red ex,varlis) else not freeofl(mvar ex,varlis) or nonlnr(lc ex,varlis) or nonlnr(red ex,varlis); % ***** Support for one_of and root_of *****. symbolic procedure mkrootsoftag(); begin scalar name; integer n; loop: n:=n #+1; name := intern compress append('(t a g _),explode n); if flagp(name,'used!*) then go to loop; return reval name; end; 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 !*k2q list('root_of,e1,x,name),list var,mu) end; put('root_of,'psopfn,'root_of_eval); symbolic procedure root_of_eval u; begin scalar !*!*norootvarrenamep!*!*,x,n; n := if cddr u then caddr u else mkrootsoftag(); !*!*norootvarrenamep!*!* := n; x := solveeval{car u,cadr u}; if eqcar(x,'list) then x := cdr x else typerr(x,"list"); x := foreach j in x collect if eqcar(j,'equal) then caddr j else typerr(j,"equation"); if null x then rederr "solve confusion in root_of_eval" else if null cdr x then return car x else return{'one_of, 'list . x,n} end; put('root_of,'subfunc,'subrootof); symbolic procedure subrootof(l,expn); % Sets up a formal SUB expression when necessary. begin scalar x,y; for each j in cddr expn do if (x := assoc(j,l)) then <<y := x . y; l := delete(x,l)>>; expn := sublis(l,car expn) . for each j in cdr expn collect subsublis(l,j); %to ensure only opr and individual args are transformed; if null y then return expn; expn := aconc!*(for each j in reversip!* y collect list('equal,car j,aeval cdr j),expn); return if l then subeval expn else mk!*sq !*p2q mksp('sub . expn,1) end; %-(algebraic << %- %-depend(!~p,!~x); % Needed for the simplification of the rule pattern. %- %-let root_of(~p,~x,~tg)^~n => %- sub(x=root_of(p,x,tg), %- -reduct(p,x)/coeffn(p,x,deg(p,x))) ^ (n-deg(p,x)+1) %- when fixp n and deg(p,x)>=1 and n>=deg(p,x); %- %-nodepend(!~p,!~x); %- %->>) where dmode!*=nil,!*modular=nil,!*rounded=nil,!*complex=nil; symbolic procedure polyp(p,f)$ % tests if p is a polynomial in f and its derivatives % p: expression % f: function if my_freeof(p,f) then t else begin scalar a$ if atom p then a:=t else if member(car p,list('EXPT,'PLUS,'MINUS,'TIMES,'QUOTIENT,'DF)) then % erlaubte Funktionen <<if (car p='PLUS) or (car p='TIMES) then <<p:=cdr p$ while p do if a:=polyp(car p,f) then p:=cdr p else p:=nil>> else if (car p='MINUS) then a:=polyp(cadr p,f) else if (car p='QUOTIENT) then <<if freeof(caddr p,f) then a:=polyp(cadr p,f)>> else if car p='EXPT then % Exponent <<if (fixp caddr p) then if caddr p>0 then a:=polyp(cadr p,f)>> else if car p='DF then % Ableitung if (cadr p=f) or freeof(cadr p,f) then a:=t>> else a:=(p=f)$ return a end$ symbolic procedure polypeval u; begin scalar bool,v; v := cadr u; u := simpcar u; if cdr u neq 1 then return nil else u := kernels car u; while u and null bool do <<if v neq car u and smember(v,car u) then bool := t; u := cdr u>>; return null bool end; put('polyp,'psopfn,'polypeval); (algebraic << depend(!~p,!~x); let root_of(~p,~x,~tg)^~n => sub(x=root_of(p,x,tg), -reduct(p,x)/coeffn(p,x,deg(p,x))) ^ (n-deg(p,x)+1) when polyp(p,x) and fixp n and deg(p,x)>=1 and n>=deg(p,x); nodepend(!~p,!~x); >>) where dmode!*=nil,!*modular=nil,!*rounded=nil,!*complex=nil; symbolic procedure expand_cases u; begin scalar bool,sl,tags; sl:=list nil; tags:=list nil; u := reval u; if not eqcar(u,'list) then typerr(u,"equation list") else u := cdr u; if eqcar(car u,'list) then <<u := for each j in u collect if eqcar(j,'list) then cdr j else typerr(j,"equation list"); bool := t>> else u := for each j in u collect {j}; u := for each j in u join expand_case1(j,sl,tags); return 'list . for each j in u collect if null bool then car j else 'list . j end; symbolic procedure expand_case1(u,sl,tags); if null u then nil else expand_merge(expand_case2(car u,sl,tags), expand_case1(cdr u,sl,tags)); symbolic procedure expand_merge(u,v); if null v then for each j in u collect {j} else for each j in u join for each k in v collect j . k; symbolic procedure expand_case2(u,sl,tags); begin scalar tag,v,var; var := cadr u; v := caddr u; if eqcar(v,'one_of) then <<tag := caddr v; if tag member tags then typerr(tag,"unique choice tag") else if null assoc(tag,sl) then cdr sl := (tag . cdadr v) . cdr sl; return if eqcar(cadr v,'list) then for each j in cdadr v collect {'equal,var,j} else typerr(cadr v,"list")>> % The next section doesn't do anything currently since root_of % is wrapped in a !*sq at this point. else if eqcar(v,'root_of) then <<tag := cadddr v; cdr tags := tag . cdr tags; if assoc(tag,sl) then typerr(tag,"unique choice tag")>>; return {u} end; % Rules for solving inverse trigonometrical functions. fluid '(solve_invtrig_soln!*); share solve_invtrig_soln!*; symbolic procedure check_solve_inv_trig(fn,equ,var); begin scalar x,s; x := evalletsub2({'(solve_trig_rules),{'simp!*,mkquote {fn,equ}}}, nil); if errorp x or not ((x := car x) freeof '(asin acos atan)) then return nil; for each sol in cdr solveeval {mk!*sq subtrsq(x,simp!* {fn,0}), var} do if is_solution(sol,equ) then s := caddr sol . s; if null s then <<solve_invtrig_soln!* := 1; return t>> % no solution found else if null cdr s then s := car s % one solution else s := 'one_of . s; solve_invtrig_soln!* := {'difference,var,s}; return t end; flag('(check_solve_inv_trig),'boolean); symbolic procedure is_solution(sol,equ); begin scalar var,s,rhs,result; var := cadr sol; rhs := caddr sol; equ := numr simp!* equ; if eqcar(rhs,'one_of) then result := check!-solns(for each s in cdr rhs collect {{simp!* s},{var},1}, equ,var) else if eqcar(rhs,'root_of) then result := t else result := check!-solns({{{simp!* rhs},{var},1}},equ,var); return if not (result eq 'unsolved) then result else nil end; symbolic procedure check!-condition u; null !*precise or eval formbool(u,nil,'algebraic); endmodule; module ppsoln; % Solve surd eqns, mainly by principle of powers method. % Authors: Anthony C. Hearn and Stanley L. Kameny. fluid '(!*complex !*msg !*numval !*ppsoln); global '(bfone!*); !*ppsoln := t; % Keep this as internal switch. 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!-solns(z,numr w,var); % return if z eq 'unsolved then list list(list w,nil,mu) else z return if z eq 'unsolved then mkrootsof(w,var,mu) else z end; symbolic procedure principle!-of!-powers!-soln(ex,x1,var,mu); % Finds solutions of ex=0 by the principle of powers method. Return % 'unsolved if solutions can't be found. begin scalar z; 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 := updkorder 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!-solns(z,ex,var); return z end; symbolic procedure check!-solns(z,ex,var); begin scalar x,y,fv,sx,vs; fv := freevarl(ex,var); for each z1 in z do fv := union(fv,union(freevarl(numr caar z1,var), freevarl(denr caar z1,var))); fv := delete('i,fv); % this does only one random setting!! if fv then for each v in fv do if not flagp(v,'constant) then vs := (v . list('quotient,1+random 999,1000)) . vs; sx := if vs then numr subf(ex,vs) else ex; while z do if null cadar z then <<z := nil; x := 'unsolved>> else if <<y := numr subf(ex,list(caadar z . mk!*sq caaar z)); null y % to do multiple random tests, the vs, sx setting and testing % would be moved here and done in a loop. or fv and null(y := numr subf(sx,list(caadar z . mk!*sq subsq(caaar z,vs)))) or null numvalue y>> then <<x := car z . x; z := cdr z>> else z := cdr z; return if null x then 'unsolved else x end; symbolic procedure suitable!-expt u; eqcar(u,'expt) and eqcar(caddr u,'quotient) and cadr caddr u = 1 and fixp caddr caddr u; symbolic procedure freevarl(ex,var); <<for each k in allkern list(ex ./ 1) do l := union(l,varsift(k,var)); delete(var,l)>> where l=if var then list var else nil; 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) 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); symbolic procedure numvalue u; % Find floating point value of sf u. begin scalar !*numval,x,c,cp,p,m; m := !*msg; !*msg := nil; !*numval := t; c := ('i memq freevarl(u,nil)); if (cp := !*complex) then off complex; x := setdmode('rounded,t); p := precision 10; if x eq '!:rd!: then x := 'rounded; % <==== to avoid error later if c then on complex; !*msg := m; u := numr simp prepf u; !*msg := nil; if c then off complex; if x then setdmode(x,t) else setdmode('rounded,nil); if cp then on complex; precision p; !*msg := m; return if eqcar(u,'!:rd!:) and (numvchk(100,z) where z=round!* u) or eqcar(u,'!:cr!:) and (numvchk(10,z) where z=retag crrl u) and (numvchk(10,z) where z=retag crim u) then nil else u end; symbolic procedure numvchk(fact,z); if atom z then fact*abs z<1 else lessp!:(timbf(bfloat fact,abs!: z),bfone!*); endmodule; module solvelnr; % Code for solving a general system of linear eqns. % Authors: Anthony C. Hearn and Eberhard Schruefer. % Modifications by: David Hartley. % Based on code by David R. Stoutemyer modified by Donald R. Morrison. % Copyright (c) 1993 RAND. All rights reserved. % The number of equations and the number of unknowns are arbitrary. % I.e. the system can be under- or overdetermined. fluid '(!*cramer !*exp !*solvesingular asymplis!* wtl!* !*arbvars !*trsparse !*varopt bareiss!-step!-size!*); % switch solveinconsistent; % !*solveinconsistent := t; % Default value. symbolic procedure solvelnrsys(exlis,varlis); % exlis: list of sf, varlis: list of kernel % -> solvelnrsys: tagged solution list % Check the system for sparsity, then decide whether to use the % Cramer or Bareiss method. Using the Bareiss method on sparse % systems, 4-step elimination seems to be faster than 2-step. % The Bareiss code is not good at handling surds at the moment, % hence exptexpflistp test. begin scalar w,method; if w := solvesparsecheck(exlis,varlis) then exlis := w else exlis := exlis . varlis; if null !*cramer and null exptexpflistp 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; symbolic procedure exptexpflistp u; % True if any of u contains an expt kernel. u and (exptexpfp car u or exptexpflistp cdr u); symbolic procedure solvesyspost(exlis,varlis); % exlis: tagged solution list, varlis: list of kernel % -> solvesyspost: tagged solution list % Insert arbitrary constants and present % solutions in same order as in varlis. % Also reorders expressions to prevailing kernel order. car exlis . foreach s in cdr exlis collect if car s and null cadr s then s else begin scalar arbvars,z; if !*arbvars or (null cadr s and length varlis = 1) then arbvars := foreach v in setdiff(varlis,cadr s) collect v . mvar makearbcomplex() else varlis := intersection(varlis,cadr s); z := pair(cadr s,sublis(arbvars,car s)); z := append(z,foreach p in arbvars collect car p . !*k2q cdr p); return {foreach v in varlis collect reordsq cdr atsoc(v,z), varlis,caddr s}; end; symbolic procedure solvecramer(exlis,varlis); % exlis: list of sf, varlis: list of kernel % -> solvecramer: tagged solution list % Just a different name at the moment. glnrsolve(exlis,varlis); symbolic procedure solvesparsecheck(sys,vl); % sys: list of sf, vl: list of kernel % -> solvesparsecheck: nil or {list of sf,list of kernel} % This program checks for a sparse linear system. If the % system is sparse enough, it returns (exlis.varlis) reordered % such that a maximum triangular upper diagonal form is % established. Otherwise the result is NIL. begin scalar vl1,xl,sys1,q,x,y; integer sp; % First collect a list vl1 where each variable is followed % by the number of equations where it occurs, and then % by the number of other variables which occur in these % equations (connectivity). At the same time, collect a measure % of the sparsity. sp:=0; vl1:= for each x in vl collect x . 0 . nil; foreach q in sys do foreach x in (xl := intersection(topkerns q,vl)) do <<y := assoc(x,vl1); cdr y := (cadr y + 1) . union(xl,cddr y); sp := sp + 1>>; foreach p in vl1 do cddr p := length cddr p - 1; % could drop the -1 % Drop out if density > 80% if sp > length sys * length vl * 0.8 then <<if !*trsparse then prin2t "System is not very sparse"; return nil>>; % Sort variables first by least occurrences and then least % connectivity. % Reset kernel order and reorder equations. if !*trsparse then solvesparseprint("Original sparse system",sys,vl); vl1:=foreach x in if not !*varopt then vl1 else sort(vl1,function (lambda(x,y); cadr x<cadr y or cadr x=cadr y and cddr x < cddr y)) collect car x; if !*varopt then <<foreach k in reverse vl1 do updkorder k; sys := for each q in sys collect reorder q>>; % Next sort equations in ascending order of their first variable % and then descending order of the next variable. sys1:= (nil . nil) . foreach x in vl1 collect x . nil; foreach q in sys do <<if domainp q or not member(mvar q,vl1) then y := assoc(nil,sys1) else y := assoc(mvar q,sys1); cdr y := q . cdr y>>; foreach p in cdr sys1 do if cdr p then cdr p := sort(cdr p, function solvesparsesort); % Finally split off a leading diagonal system and push the remaining % equations down. sys := nconc(foreach p in sys1 join if cdr p then {cadr p}, reversip foreach p in sys1 join if cdr p then cddr p); if !*trsparse then solvesparseprint("Variables and equations rearranged",sys,vl1); return sys.vl1; end; symbolic procedure solvesparseprint(text,sys,vl); <<terpri(); prin2t text; for each e in sys do << e := topkerns e; for each x in vl do if memq(x,e) then prin2 "*" else prin2 "-"; terpri()>>>>; symbolic procedure topkerns u; % u:sf -> topkerns:list of kernel % kernels in top level of u if domainp u then nil else mvar u . topkerns red u; symbolic procedure solvesparsesort(x,y); % x,y: sf -> solvesparsesort: bool if domainp x then nil else if domainp y then t else if mvar x = mvar y then solvesparsesort(red y,red x) else if ordop(mvar x,mvar y) then t else nil; endmodule; module glsolve; % Routines for solving a general system of linear % equations using Cramer's rule, realized through % exterior multiplication. % Author: Eberhard Schruefer. % Modifications by: D. Hartley and R.W. Tucker. % The number of equations and the number of unknowns are arbitrary. % I.e. the system can be under- or overdetermined. fluid '(!*solvesingular vars!*); % !*solveinconsistent global '(!!arbint assumptions requirements); symbolic procedure glnrsolve(u,v); % glnrsolve(u: list of sf's, v: list of kernels) % -> (xprs: list of sq's, flg: boolean) % Adapted by D. Hartley and R.W. Tucker from E. Schruefer's glnrsolve. % The equations u must be ordered with respect to the kernels v begin scalar sgn,x,y,cnds; if null u then go to b; a: x := !*sf2ex(car u,v); if null x then u := cdr u else if inconsistency!-chk x then <<cnds := car u . cnds; x := nil; u := cdr u>>; if u and null x then go to a; b: if null u then % no consistent non-zero equations if cnds then go to d else return t . {{nil,nil,1}}; % all equations were zero if null(u := cdr u) then go to d; c: if y := subs2chkex extmult(!*sf2ex(car u,v),x) then if inconsistency!-chk y then cnds := numr cancel(lc y ./ lc x) . cnds else <<assumptions := 'list . mk!*sq !*f2q lc y . (pairp assumptions and cdr assumptions); x := y>>; if (u := cdr u) then go to c; d: for each j in cnds do requirements := 'list . mk!*sq !*f2q j . (pairp requirements and cdr requirements); if cnds then return 'inconsistent . nil; if setdiff(v,lpow x) and not !*solvesingular then return 'singular . {}; if null red x then return t . {{for each j in lpow x collect nil ./ 1,lpow x,1}}; y := lc x; sgn := evenp length lpow x; u := foreach j in lpow x collect (if (sgn := not sgn) then negf f else f) where f = !*ex2sf innprodpex(delete(j,lpow x),red x); return t . {{foreach f in u collect cancel(f ./ y),lpow x,1}}; end; symbolic procedure inconsistency!-chk u; null u or ((nil memq lpow u) and inconsistency!-chk red u); endmodule; 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. 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 ); 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:='(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 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 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 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 := solveeval list(list ('equal,s,'!#), list('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; module solvetab; % Simplification rules for SOLVE. % Author: David R. Stoutemyer. % Modifications by: Anthony C. Hearn, Donald R. Morrison, Rainer % Schoepf and Herbert Melenk. put('asin, 'inverse, 'sin); put('acos, 'inverse, 'cos); put('atan,'inverse,'tan); put('acot,'inverse,'cot); put('asec,'inverse,'sec); put('acsc,'inverse,'csc); algebraic; Comment Rules for reducing the number of distinct kernels in an equation; operator sol; % for all a,b,c,d,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 not fixp c and ratnump c and not fixp d and ratnump d let sol(a**c-b**d, x) = a**(c*lcm(den c,den d)) - b**(d*lcm(den c,den 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) = if !*allbranch then sin((a+b)/2)* cos((a-b)/2) else a+b, 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) = cos((a+b)/2)*cos((a-b)/2), 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; solve_trig_rules := {sin(~x + ~y) => sin x * cos y + cos x * sin y, sin(~x - ~y) => sin x * cos y - cos x * sin y, cos(~x + ~y) => cos x * cos y - sin x * sin y, cos(~x - ~y) => cos x * cos y + sin x * sin y}; fluid '(solve_invtrig_soln!*); share solve_invtrig_soln!*; clear solve_invtrig_soln!*; invtrig_solve_rules := { sol(asin(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('sin,asin(x) + y,z), sol(acos(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('cos,acos(x) + y,z), sol(atan(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('tan,atan(x) + y,z), sol(acos(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('sin,acos(x) + y,z), sol(atan(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('sin,atan(x) + y,z), sol(asin(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('cos,asin(x) + y,z), sol(atan(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('cos,atan(x) + y,z), sol(~n*asin(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('sin,n*asin(x) + y,z), sol(~n*acos(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('cos,n*acos(x) + y,z), sol(~n*acos(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('sin,n*acos(x) + y,z), sol(~n*atan(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('sin,n*atan(x) + y,z), sol(~n*asin(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('cos,n*asin(x) + y,z), sol(~n*atan(~x) + ~y,~z) => solve_invtrig_soln!* when check_solve_inv_trig('cos,n*atan(x) + y,z) }; let invtrig_solve_rules; % The following rules allow REDUCE to solve some classes of equations % where a variable appears inside and outside a log or an exponential. % The results are based on Lambert's W (Omega) function which is fully % supported in the specfn package. The ruleset has one central rule % which produces the Omega function expression in the simplest (rather % special) form, while the more general cases are mapped towards this % rule by reforming the equation algebraically or by variable % transformations. lambert_rules := { % Basic solution of x=log(c*x/d) sol(~x + log(~~c*~x/~~d),~x) => x - lambert_w(d/c) when c freeof x and d freeof x, % General forms transformed to simpler ones. sol(~~a*~x + ~~b*log(~c) + ~w,x) => sol(a*x + b*log(c*e^(w/b)), x) when a freeof x and b freeof x and w freeof x and not(c freeof x), sol(~~a*~x + ~~b*log(~~c*x/~~d),x) => sub(x=a*x/b, sol(x + log(c*b*x/(a*d)),x)) when (a neq 1 or b neq 1) and a freeof x and b freeof x and c freeof x and d freeof x, sol(~~a*~x + ~~b*log((~~c*x + ~u)/~~d),x) => sub(x=x+u/c, sol(num(a*(x-u/c) + b*log(c*x/d)),x)) when a freeof x and b freeof x and c freeof x and d freeof x and u freeof x, sol(~~a*~x + ~~b*log((~~c*x^~n)/~~d),x) => sol(num(a*x + n*b*log x + 1/n*log(c/d)),x) when a freeof x and b freeof x and c freeof x and d freeof x and n freeof x, sol(~~a*~x^~~n + ~~b*e^(~~c*~x/~~d),x) => sol(num(log(a) + n*log(x) - (log(-b)*d + c*x)/d), x) when a freeof x and b freeof x and c freeof x and d freeof x and n freeof x, sol(~~a*~x + ~~b*e^(~~c*~x/~~d) + ~f,x) => sub(x=a*x+f/a,sol(num(x + b*e^(-c*f/(a*d))*e^(c*x/(a*d))),x)) when a freeof x and b freeof x and c freeof x and d freeof x and f freeof x }$ % let lambert_rules; symbolic procedure lambertp(e1,x); smemq('log,e1) or smemq('expt,e1); symbolic; fluid '(sol!-rulesets!*); sol!-rulesets!*:={{'lambertp,'lambert_rules}}; symbolic procedure solve!-apply!-rules(e1,var); begin scalar rules,u; u:=list('sol,mk!*sq(e1 ./ 1), var); for each r in sol!-rulesets!* do if apply(car r,{e1,var}) then rules := cadr r . rules; if null rules then return simp!* u; return car evalletsub2({rules,{'simp!*, mkquote u}},nil); end; endmodule; module quartic; % Procedures for solving cubic, quadratic and quartic % eqns. % Author: Anthony C. Hearn. % Modifications by: Stanley L. Kameny. fluid '(!*sub2 !*rounded !*trigform dmode!*); !*trigform := t; % Default value. switch trigform; 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; % Rationalizing the value in this and the following function leads % usually to neater results. % rationalizesq simpexpt list(mk!*sq subs2!* u,'(quotient 1 3)); % simprad(u,3); symbolic procedure sqrtq u; % rationalizesq 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. % Modified to use root_val to compute numeric roots. SLK. if !*rounded and numcoef a0 and numcoef a1 and numcoef a2 then for each z in cdr root_val list mkpolyexp2(a2,a1,a0) collect simp!* z else 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 numcoef a; denr a = 1 and domainp numr a; symbolic procedure mkpolyexp2(a2,a1,a0); % The use of 'x is arbitrary here, since it is not used by root_val. <<a0 := numr a0; if numr a1 then a0 := (('x . 1) . numr a1) . a0; mk!*sq(((('x . 2) . numr a2) . a0) . 1)>>; 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)); if null numr q or not !*trigform or not all_real(a0,a1,a2) then go to cbr; % this section uses trig functions, but only when a0,a1,a2 are real. x := sqrtq negsq addsq(exptsq(q,3),exptsq(r,2)); if one_real simp list('times,'i,mk!*sq x) and not pos_num q then x := negsq x; s1 := quotsqf(atan2q(x,r),3); s2 := negsq sqrtq negsq q; sp := negsq multfq(2,multsq(s2,cossq s1)); sm := multsq(simp '(sqrt 3),multsq(s2,sinsq s1)); % sp := -2*sqrt(-q)*cos(atan2(sqrt( - q**3 - r**2),r)/3)$ % sm := - sqrt(-q)*sqrt(3)*sin(atan2(sqrt( - q**3 - r**2),r)/3)$ go to com; cbr: 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); com: 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 pos_num a; begin scalar r,dmode,!*msg,!*numval; dmode := dmode!*; !*numval := t; on rounded,complex; a := resimp a; a := real_1 a and (numr simp list('sign,mk!*sq a)=1); off rounded,complex; if dmode then onoff(get(dmode,'dname),t); return a end; symbolic procedure sinsq a; simpiden list('sin,mk!*sq subs2!* a); symbolic procedure cossq a; simpiden list('cos,mk!*sq subs2!* a); symbolic procedure all_real(a,b,c); begin scalar r,dmode,!*msg,!*numval; dmode := dmode!*; !*numval := t; on complex,rounded; a := resimp a; b := resimp b; c := resimp c; a := real_1 a and real_1 b and real_1 c; off rounded,complex; if dmode then onoff(get(dmode,'dname),t); return a end; symbolic procedure real_1 x; numberp denr x and domainp numr x and null numr impartsq x; symbolic procedure one_real a; begin scalar r,dmode,!*msg,!*numval; dmode := dmode!*; !*numval := t; on complex,rounded; a := real_1 resimp a; off rounded,complex; if dmode then onoff(get(dmode,'dname),t); return a end; symbolic procedure atan2q(b,a); % Used by solvecubic to set up trig form expressions for atan2 in % terms of atan and, where necessary, a bias of +/- pi; or to call % atan2 directly if numerical solution is called for. ((begin scalar !*msg,x,y,r,dmode,q,fg,s1,s2,s3,s4,s5; y := b := simp!*(b := mk!*sq subs2!* b); x := a := simp!*(a := mk!*sq subs2!* a); if domainp numr y and domainp numr x and numberp denr y and numberp denr x then go to aret; dmode := dmode!*; on complex,rounded; y := resimp b; x := resimp a; if not(domainp numr y and domainp numr x and numberp denr y and numberp denr x) then go to ret; q := sqrtq addsq(exptsq(x,2),exptsq(y,2)); x := quotsq(x,q); y := quotsq(y,q); if null numr x then <<s1 := t; if numr simp list('sqn,list('repart,mk!*sq y))=-1 then s2 := t; go to ret>>; s3 := t; x := numr simp list('sign,list('repart,mk!*sq x)); if x=-1 then <<y := numr simp list('sign,list('repart,mk!*sq y)); if y=-1 then s4 := t else s5 := t>>; ret: off rounded,complex; if dmode then onoff(get(dmode,'dname),t); if s1 then fg := quotsqf(simp 'pi,2); if s2 then fg := negsq fg; if s3 then fg := simpiden list('atan,mk!*sq quotsq(b,a)); if s4 then fg := subtrsq(fg,simp 'pi); if s5 then fg := addsq(fg,simp 'pi); aret: return if fg then fg else simpiden list('atan2,mk!*sq b,mk!*sq a) end) where !*numval=t); 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,yy,cx,z,s,l,zz1,zz2,r,dmode,neg,!*msg,!*numval; % Convert equation to monomial form. dmode := dmode!*; 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 the % real root if there is only one; or if there are three, we choose % one that yields real coefficients for the quadratics. If no % roots are known to be real, we use an arbitrary one. yy := subtrsq(exptsq(a3,2),multfq(4,a2)); x := solvecubic(!*f2q 1, negsq a2, subs2!* subtrsq(multsq(a1,a3),multfq(4,a0)), subs2!* negsq addsq(exptsq(a1,2), multsq(a0,yy))); cx := car x; % Now check for real roots of the cubic. for each rr in x do if one_real rr then s := append(s,list rr); x := if (l := length s)=1 then car s else cx; % Now solve the two equivalent quadratic equations. a3 := quotsqf(a3,2); yy := quotsqf(yy,4); % select real coefficient for quadratic if possible. y := addsq(yy,x); if l<2 then go to zz; loop: if not pos_num negsq y then go to zz else if l=1 then <<x := cx; y := addsq(yy,x); go to zz>>; l := l-1; s := cdr s; x := car s; y := addsq(yy,x); go to loop; zz: y := sqrtq y; x := quotsqf(x,2); z := sqrtq subtrsq(exptsq(x,2),a0); % the following test is needed, according to some editions of % Abramowitz and Stegun, to select the correct signs % (for the terms z) in the quadratics to produce correct roots. % Unfortunately, this test may fail for coefficients which are not % numeric because of the inability to recognize zero. !*numval := t; on rounded,complex; if null numr (zz1 := resimp subtrsq(a1,addsq(multsq(subtrsq(a3,y),addsq(x,z)), multsq(addsq(a3,y),subtrsq(x,z))))) then go to rst; if null numr (zz2 := resimp subtrsq(a1,addsq(multsq(subtrsq(a3,y),subtrsq(x,z)), multsq(addsq(a3,y),addsq(x,z))))) then <<neg := t; go to rst>>; if domainp numr zz1 and domainp numr zz2 and numberp denr zz1 and numberp denr zz2 and numr simp list('sign,list('difference,list('norm,mk!*sq zz1), list('norm,mk!*sq zz2)))=1 then neg := t; rst: off rounded,complex; if dmode then onoff(get(dmode,'dname),t); if neg then z := negsq z; return append(solvequadratic(!*f2q 1,subtrsq(a3,y),subtrsq(x,z)), solvequadratic(!*f2q 1,addsq(a3,y),addsq(x,z))) end; endmodule; end;