\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra solvefor.spad} \author{Stephen M. Watt, Barry Trager} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package SOLVEFOR PolynomialSolveByFormulas} <>= )abbrev package SOLVEFOR PolynomialSolveByFormulas -- Current fields with **: (%, RationalNumber) -> % are -- ComplexFloat, RadicalExtension(K) and RationalRadical -- SMW June 86, BMT Sept 93 ++ Description: ++ This package factors the formulas out of the general solve code, ++ allowing their recursive use over different domains. ++ Care is taken to introduce few radicals so that radical extension ++ domains can more easily simplify the results. PolynomialSolveByFormulas(UP, F): PSFcat == PSFdef where UP: UnivariatePolynomialCategory F F: Field with **: (%, Fraction Integer) -> % L ==> List PSFcat == with solve: UP -> L F ++ solve(u) \undocumented particularSolution: UP -> F ++ particularSolution(u) \undocumented mapSolve: (UP, F -> F) -> Record(solns: L F, maps: L Record(arg:F,res:F)) ++ mapSolve(u,f) \undocumented linear: UP -> L F ++ linear(u) \undocumented quadratic: UP -> L F ++ quadratic(u) \undocumented cubic: UP -> L F ++ cubic(u) \undocumented quartic: UP -> L F ++ quartic(u) \undocumented -- Arguments give coefs from high to low degree. linear: (F, F) -> L F ++ linear(f,g) \undocumented quadratic: (F, F, F) -> L F ++ quadratic(f,g,h) \undocumented cubic: (F, F, F, F) -> L F ++ cubic(f,g,h,i) \undocumented quartic: (F, F, F, F, F) -> L F ++ quartic(f,g,h,i,j) \undocumented aLinear: (F, F) -> F ++ aLinear(f,g) \undocumented aQuadratic: (F, F, F) -> F ++ aQuadratic(f,g,h) \undocumented aCubic: (F, F, F, F) -> F ++ aCubic(f,g,h,j) \undocumented aQuartic: (F, F, F, F, F) -> F ++ aQuartic(f,g,h,i,k) \undocumented PSFdef == add ----------------------------------------------------------------- -- Stuff for mapSolve ----------------------------------------------------------------- id ==> (IDENTITY$Lisp) maplist: List Record(arg: F, res: F) := [] mapSolving?: Boolean := false -- map: F -> F := id #1 replaced with line below map: Boolean := false mapSolve(p, fn) == -- map := fn #1 replaced with line below locmap: F -> F := fn #1; map := id locmap mapSolving? := true; maplist := [] slist := solve p mapSolving? := false; -- map := id #1 replaced with line below locmap := id #1; map := id locmap [slist, maplist] part(s: F): F == not mapSolving? => s -- t := map s replaced with line below t: F := SPADCALL(s, map)$Lisp t = s => s maplist := cons([t, s], maplist) t ----------------------------------------------------------------- -- Entry points and error handling ----------------------------------------------------------------- cc ==> coefficient -- local intsolve intsolve(u:UP):L(F) == u := (factors squareFree u).1.factor n := degree u n=1 => linear (cc(u,1), cc(u,0)) n=2 => quadratic (cc(u,2), cc(u,1), cc(u,0)) n=3 => cubic (cc(u,3), cc(u,2), cc(u,1), cc(u,0)) n=4 => quartic (cc(u,4), cc(u,3), cc(u,2), cc(u,1), cc(u,0)) error "All sqfr factors of polynomial must be of degree < 5" solve u == ls := nil$L(F) for f in factors squareFree u repeat lsf := intsolve f.factor for i in 1..(f.exponent) repeat ls := [:lsf,:ls] ls particularSolution u == u := (factors squareFree u).1.factor n := degree u n=1 => aLinear (cc(u,1), cc(u,0)) n=2 => aQuadratic (cc(u,2), cc(u,1), cc(u,0)) n=3 => aCubic (cc(u,3), cc(u,2), cc(u,1), cc(u,0)) n=4 => aQuartic (cc(u,4), cc(u,3), cc(u,2), cc(u,1), cc(u,0)) error "All sqfr factors of polynomial must be of degree < 5" needDegree(n: Integer, u: UP): Boolean == degree u = n => true error concat("Polynomial must be of degree ", n::String) needLcoef(cn: F): Boolean == cn ~= 0 => true error "Leading coefficient must not be 0." needChar0(): Boolean == characteristic$F = 0 => true error "Formula defined only for fields of characteristic 0." linear u == needDegree(1, u) linear (coefficient(u,1), coefficient(u,0)) quadratic u == needDegree(2, u) quadratic (coefficient(u,2), coefficient(u,1), coefficient(u,0)) cubic u == needDegree(3, u) cubic (coefficient(u,3), coefficient(u,2), coefficient(u,1), coefficient(u,0)) quartic u == needDegree(4, u) quartic (coefficient(u,4),coefficient(u,3), coefficient(u,2),coefficient(u,1),coefficient(u,0)) ----------------------------------------------------------------- -- The formulas ----------------------------------------------------------------- -- local function for testing equality of radicals. -- This function is necessary to detect at least some of the -- situations like sqrt(9)-3 = 0 --> false. equ(x:F,y:F):Boolean == ( (recip(x-y)) case "failed" ) => true false linear(c1, c0) == needLcoef c1 [- c0/c1 ] aLinear(c1, c0) == first linear(c1,c0) quadratic(c2, c1, c0) == needLcoef c2; needChar0() (c0 = 0) => [0$F,:linear(c2, c1)] (c1 = 0) => [(-c0/c2)**(1/2),-(-c0/c2)**(1/2)] D := part(c1**2 - 4*c2*c0)**(1/2) [(-c1+D)/(2*c2), (-c1-D)/(2*c2)] aQuadratic(c2, c1, c0) == needLcoef c2; needChar0() (c0 = 0) => 0$F (c1 = 0) => (-c0/c2)**(1/2) D := part(c1**2 - 4*c2*c0)**(1/2) (-c1+D)/(2*c2) w3: F := (-1 + (-3::F)**(1/2)) / 2::F cubic(c3, c2, c1, c0) == needLcoef c3; needChar0() -- case one root = 0, not necessary but keeps result small (c0 = 0) => [0$F,:quadratic(c3, c2, c1)] a1 := c2/c3; a2 := c1/c3; a3 := c0/c3 -- case x**3-a3 = 0, not necessary but keeps result small (a1 = 0 and a2 = 0) => [ u*(-a3)**(1/3) for u in [1, w3, w3**2 ] ] -- case x**3 + a1*x**2 + a1**2*x/3 + a3 = 0, the general for- -- mula is not valid in this case, but solution is easy. P := part(-a1/3::F) equ(a1**2,3*a2) => S := part((- a3 + (a1**3)/27::F)**(1/3)) [ P + S*u for u in [1,w3,w3**2] ] -- general case Q := part((3*a2 - a1**2)/9::F) R := part((9*a1*a2 - 27*a3 - 2*a1**3)/54::F) D := part(Q**3 + R**2)**(1/2) S := part(R + D)**(1/3) -- S = 0 is done in the previous case [ P + S*u - Q/(S*u) for u in [1, w3, w3**2] ] aCubic(c3, c2, c1, c0) == needLcoef c3; needChar0() (c0 = 0) => 0$F a1 := c2/c3; a2 := c1/c3; a3 := c0/c3 (a1 = 0 and a2 = 0) => (-a3)**(1/3) P := part(-a1/3::F) equ(a1**2,3*a2) => S := part((- a3 + (a1**3)/27::F)**(1/3)) P + S Q := part((3*a2 - a1**2)/9::F) R := part((9*a1*a2 - 27*a3 - 2*a1**3)/54::F) D := part(Q**3 + R**2)**(1/2) S := part(R + D)**(1/3) P + S - Q/S quartic(c4, c3, c2, c1, c0) == needLcoef c4; needChar0() -- case one root = 0, not necessary but keeps result small (c0 = 0) => [0$F,:cubic(c4, c3, c2, c1)] -- Make monic: a1 := c3/c4; a2 := c2/c4; a3 := c1/c4; a4 := c0/c4 -- case x**4 + a4 = 0 <=> (x**2-sqrt(-a4))*(x**2+sqrt(-a4)) -- not necessary but keeps result small. (a1 = 0 and a2 = 0 and a3 = 0) => append( quadratic(1, 0, (-a4)**(1/2)),_ quadratic(1 ,0, -((-a4)**(1/2))) ) -- Translate w = x+a1/4 to eliminate a1: w**4+p*w**2+q*w+r p := part(a2-3*a1*a1/8::F) q := part(a3-a1*a2/2::F + a1**3/8::F) r := part(a4-a1*a3/4::F + a1**2*a2/16::F - 3*a1**4/256::F) -- t0 := the cubic resolvent of x**3-p*x**2-4*r*x+4*p*r-q**2 -- The roots of the translated polynomial are those of -- two quadratics. (What about rt=0 ?) -- rt=0 can be avoided by picking a root ~= p of the cubic -- polynomial above. This is always possible provided that -- the input is squarefree. In this case the two other roots -- are +(-) 2*r**(1/2). if equ(q,0) -- this means p is a root then t0 := part(2*(r**(1/2))) else t0 := aCubic(1, -p, -4*r, 4*p*r - q**2) rt := part(t0 - p)**(1/2) slist := append( quadratic( 1, rt, (-q/rt + t0)/2::F ), quadratic( 1, -rt, ( q/rt + t0)/2::F )) -- Translate back: [s - a1/4::F for s in slist] aQuartic(c4, c3, c2, c1, c0) == needLcoef c4; needChar0() (c0 = 0) => 0$F a1 := c3/c4; a2 := c2/c4; a3 := c1/c4; a4 := c0/c4 (a1 = 0 and a2 = 0 and a3 = 0) => (-a4)**(1/4) p := part(a2-3*a1*a1/8::F) q := part(a3-a1*a2/2::F + a1**2*a1/8::F) r := part(a4-a1*a3/4::F + a1**2*a2/16::F - 3*a1**4/256::F) if equ(q,0) then t0 := part(2*(r**(1/2))) else t0 := aCubic(1, -p, -4*r, 4*p*r - q**2) rt := part(t0 - p)**(1/2) s := aQuadratic( 1, rt, (-q/rt + t0)/2::F ) s - a1/4::F @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @ <<*>>= <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}