Initial population.

1 files changed, 327 insertions, 0 deletions

diff --git a/src/algebra/solvefor.spad.pamphlet b/src/algebra/solvefor.spad.pamphlet
new file mode 100644
index 00000000..7095d4e1
--- /dev/null
+++ b/src/algebra/solvefor.spad.pamphlet
@@ -0,0 +1,327 @@
+\documentclass{article}
+\usepackage{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}
+<<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}
+<<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.
+@
+<<*>>=
+<<license>>
+
+<<package SOLVEFOR PolynomialSolveByFormulas>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}