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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra solverad.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package SOLVERAD RadicalSolvePackage}
<<package SOLVERAD RadicalSolvePackage>>=
)abbrev package SOLVERAD RadicalSolvePackage
++ Author: P.Gianni
++ Date Created: Summer 1990
++ Date Last Updated: October 1991
++ Basic Functions:
++ Related Constructors: SystemSolvePackage, FloatingRealPackage,
++ FloatingComplexPackage
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package tries to find solutions
++ expressed in terms of radicals for systems of equations
++ of rational functions with coefficients in an integral domain R.
RadicalSolvePackage(R): Cat == Capsule where
    R   :  Join(EuclideanDomain, CharacteristicZero)
    PI ==> PositiveInteger
    NNI==> NonNegativeInteger
    Z  ==> Integer
    B  ==> Boolean
    ST ==> String
    PR ==> Polynomial R
    UP ==> SparseUnivariatePolynomial PR
    LA ==> LocalAlgebra(PR, Z, Z)
    RF ==> Fraction PR
    RE ==> Expression R
    EQ ==> Equation
    SY ==> Symbol
    SU ==> SuchThat(List RE, List Equation RE)
    SUP==> SparseUnivariatePolynomial
    L  ==> List
    P  ==> Polynomial

    SOLVEFOR ==> PolynomialSolveByFormulas(SUP RE, RE)
    UPF2     ==> SparseUnivariatePolynomialFunctions2(PR,RE)

    Cat ==> with

        radicalSolve :     (RF,SY)      -> L EQ RE
          ++ radicalSolve(rf,x) finds the solutions expressed in terms of
          ++ radicals of the equation rf = 0 with respect to the symbol x,
          ++ where rf is a rational function.
        radicalSolve :       RF         -> L EQ RE
          ++ radicalSolve(rf) finds the solutions expressed in terms of
          ++ radicals of the equation rf = 0, where rf is a
          ++ univariate rational function.
        radicalSolve :    (EQ RF,SY)    -> L EQ RE
          ++ radicalSolve(eq,x) finds the solutions expressed in terms of
          ++ radicals of the equation of rational functions eq
          ++ with respect to the symbol x.
        radicalSolve :      EQ RF       -> L EQ RE
          ++ radicalSolve(eq) finds the solutions expressed in terms of
          ++ radicals of the equation of rational functions eq
          ++ with respect to the unique symbol x appearing in eq.
        radicalSolve :    (L RF,L SY)   -> L L EQ RE
          ++ radicalSolve(lrf,lvar) finds the solutions expressed in terms of
          ++ radicals of the system of equations lrf = 0 with
          ++ respect to the list of symbols lvar,
          ++ where lrf is a list of rational functions.
        radicalSolve :       L RF       -> L L EQ RE
          ++ radicalSolve(lrf) finds the solutions expressed in terms of
          ++ radicals of the system of equations lrf = 0, where lrf is a
          ++ system of univariate rational functions.
        radicalSolve :   (L EQ RF,L SY) -> L L EQ RE
          ++ radicalSolve(leq,lvar) finds the solutions expressed in terms of
          ++ radicals of the system of equations of rational functions leq
          ++ with respect to the list of symbols lvar.
        radicalSolve :     L EQ RF      -> L L EQ RE
          ++ radicalSolve(leq) finds the solutions expressed in terms of
          ++ radicals of the system of equations of rational functions leq
          ++ with respect to the unique symbol x appearing in leq.
        radicalRoots :      (RF,SY)     -> L RE
          ++ radicalRoots(rf,x) finds the roots expressed in terms of radicals
          ++ of the rational function rf with respect to the symbol x.
        radicalRoots :    (L RF,L SY)   -> L L RE
          ++ radicalRoots(lrf,lvar) finds the roots expressed in terms of
          ++ radicals of the list of rational functions lrf
          ++ with respect to the list of symbols lvar.
        contractSolve:    (EQ RF,SY)  -> SU
          ++ contractSolve(eq,x) finds the solutions expressed in terms of
          ++ radicals of the equation of rational functions eq
          ++ with respect to the symbol x.  The result contains new
          ++ symbols for common subexpressions in order to reduce the
          ++ size of the output.
        contractSolve:    (RF,SY)     -> SU
          ++ contractSolve(rf,x) finds the solutions expressed in terms of
          ++ radicals of the equation rf = 0 with respect to the symbol x,
          ++ where rf is a rational function. The result contains  new
          ++ symbols for common subexpressions in order to reduce the
          ++ size of the output.
    Capsule ==> add
        import DegreeReductionPackage(PR, R)
        import SOLVEFOR

        SideEquations: List EQ RE := []
        ContractSoln:  B := false

        ---- Local Function Declarations ----
        solveInner:(PR, SY, B) -> SU
        linear:    UP -> List RE
        quadratic: UP -> List RE
        cubic:     UP -> List RE
        quartic:   UP -> List RE
        rad:       PI -> RE
        wrap:      RE -> RE
        New:       RE -> RE
        makeEq : (List RE,L SY) -> L EQ RE
        select :    L L RE      -> L L RE
        isGeneric? :  (L PR,L SY)  ->  Boolean
        findGenZeros :  (L PR,L SY) -> L L RE
        findZeros   :   (L PR,L SY) -> L L RE


        New s ==
            s = 0 => 0
            S := new()$Symbol ::PR::RF::RE
            SideEquations := append([S = s], SideEquations)
            S

        linear u    == [(-coefficient(u,0))::RE /(coefficient(u,1))::RE]
        quadratic u == quadratic(map(coerce,u)$UPF2)$SOLVEFOR
        cubic u     == cubic(map(coerce,u)$UPF2)$SOLVEFOR
        quartic u   == quartic(map(coerce,u)$UPF2)$SOLVEFOR
        rad n       == n::Z::RE
        wrap s      == (ContractSoln => New s; s)


        ---- Exported Functions ----


       -- find the zeros of components in "generic" position --
        findGenZeros(rlp:L PR,rlv:L SY) : L L RE ==
         pp:=rlp.first
         v:=first rlv
         rlv:=rest rlv
         res:L L RE:=[]
         res:=append([reverse cons(r,[eval(
           (-coefficient(univariate(p,vv),0)::RE)/(leadingCoefficient univariate(p,vv))::RE,
              kernel(v)@Kernel(RE),r) for vv in rlv for p in rlp.rest])
                for r in radicalRoots(pp::RF,v)],res)
         res


        findZeros(rlp:L PR,rlv:L SY) : L L RE ==
         parRes:=[radicalRoots(p::RF,v) for p in rlp for v in rlv]
         parRes:=select parRes
         res:L L RE :=[]
         res1:L RE
         for par in parRes repeat
           res1:=[par.first]
           lv1:L Kernel(RE):=[kernel rlv.first]
           rlv1:=rlv.rest
           p1:=par.rest
           while p1~=[] repeat
             res1:=cons(eval(p1.first,lv1,res1),res1)
             p1:=p1.rest
             lv1:=cons(kernel rlv1.first,lv1)
             rlv1:=rlv1.rest
           res:=cons(res1,res)
         res

        radicalSolve(pol:RF,v:SY) ==
          [equation(v::RE,r) for r in radicalRoots(pol,v)]

        radicalSolve(p:RF) ==
          zero? p =>
             error "equation is always satisfied"
          lv:=removeDuplicates
             concat(variables numer p, variables denom p)
          empty? lv => error "inconsistent equation"
          #lv>1 => error "too many variables"
          radicalSolve(p,lv.first)

        radicalSolve(eq: EQ RF) ==
          radicalSolve(lhs eq -rhs eq)

        radicalSolve(eq: EQ RF,v:SY) ==
           radicalSolve(lhs eq - rhs eq,v)

        radicalRoots(lp: L RF,lv: L SY) ==
          parRes:=triangularSystems(lp,lv)$SystemSolvePackage(R)
          parRes= list [] => []
           -- select the components in "generic" form
          rlv:=reverse lv
          rpRes:=[reverse res for res in parRes]
          listGen:= [res for res in rpRes|isGeneric?(res,rlv)]
          result:L L RE:=[]
          if listGen~=[] then
            result:="append"/[findGenZeros(res,rlv) for res in listGen]
            for res in listGen repeat
                rpRes:=delete(rpRes,position(res,rpRes))
           --  non-generic components
          rpRes = [] => result
          append("append"/[findZeros(res,rlv) for res in rpRes],
                         result)

        radicalSolve(lp:L RF,lv:L SY) ==
          [makeEq(lres,lv) for lres in radicalRoots(lp,lv)]

        radicalSolve(lp: L RF) ==
          lv:="setUnion"/[setUnion(variables numer p,variables denom p)
                          for p in lp]
          [makeEq(lres,lv) for lres in radicalRoots(lp,lv)]

        radicalSolve(le:L EQ RF,lv:L SY) ==
          lp:=[rhs p -lhs p for p in le]
          [makeEq(lres,lv) for lres in radicalRoots(lp,lv)]

        radicalSolve(le: L EQ RF) ==
          lp:=[rhs p -lhs p for p in le]
          lv:="setUnion"/[setUnion(variables numer p,variables denom p)
                          for p in lp]
          [makeEq(lres,lv) for lres in radicalRoots(lp,lv)]

        contractSolve(eq:EQ RF, v:SY)==
           solveInner(numer(lhs eq - rhs eq), v, true)

        contractSolve(pq:RF, v:SY) == solveInner(numer pq, v, true)

        radicalRoots(pq:RF, v:SY) == lhs solveInner(numer pq, v, false)


       -- test if the ideal is radical in generic position --
        isGeneric?(rlp:L PR,rlv:L SY) : Boolean ==
          "and"/[degree(f,x)=1 for f in rest rlp  for x in rest rlv]

        ---- select  the univariate factors
        select(lp:L L RE) : L L RE ==
          lp=[] => list []
          [:[cons(f,lsel) for lsel in select lp.rest] for f in lp.first]

        ---- Local Functions ----
       -- construct the equation
        makeEq(nres:L RE,lv:L SY) : L EQ RE ==
          [equation(x :: RE,r) for x in lv for r in nres]

        solveInner(pq:PR,v:SY,contractFlag:B) ==
            SideEquations := []
            ContractSoln  := contractFlag

            factors:= factors
               (factor pq)$MultivariateFactorize(SY,IndexedExponents SY,R,PR)

            constants:  List PR     := []
            unsolved:   List PR     := []
            solutions:  List RE     := []

            for f in factors repeat
                ff:=f.factor
                not member?(v, variables (ff)) =>
                    constants := cons(ff, constants)
                u := univariate(ff, v)
                t := reduce u
                u := t.pol
                n := degree u
                l: List RE :=
                    n = 1 => linear u
                    n = 2 => quadratic u
                    n = 3 => cubic u
                    n = 4 => quartic u
                    unsolved := cons(ff, unsolved)
                    []
                for s in l repeat
                    if t.deg > 1 then s := wrap s
                    T0 := expand(s, t.deg)
                    for i in 1..f.exponent repeat
                        solutions := append(T0, solutions)
                    re := SideEquations
            [solutions, SideEquations]$SU

@
\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 SOLVERAD RadicalSolvePackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}