path: root/src/algebra/solvefor.spad.pamphlet

\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}
<<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}