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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra solvelin.spad}
\author{Patrizia Gianni, Stephen M. Watt, Robert Sutor}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package LSMP LinearSystemMatrixPackage}
<<package LSMP LinearSystemMatrixPackage>>=
)abbrev package LSMP LinearSystemMatrixPackage
++ Author: P.Gianni, S.Watt
++ Date Created: Summer 1985
++ Date Last Updated:Summer 1990
++ Basic Functions: solve, particularSolution, hasSolution?, rank
++ Related Constructors: LinearSystemMatrixPackage1
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package solves linear system in the matrix form \spad{AX = B}.

LinearSystemMatrixPackage(F, Row, Col, M): Cat == Capsule where
    F: Field
    Row: Join(FiniteLinearAggregate F,ShallowlyMutableAggregate F)
    Col: Join(FiniteLinearAggregate F,ShallowlyMutableAggregate F)
    M  : MatrixCategory(F, Row, Col)

    N        ==> NonNegativeInteger
    PartialV ==> Union(Col, "failed")
    Both     ==> Record(particular: PartialV, basis: List Col)

    Cat ==> with
        solve       : (M, Col) -> Both
          ++  solve(A,B) finds a particular solution of the system \spad{AX = B}
          ++  and a basis of the associated homogeneous system \spad{AX = 0}.
        solve       : (M, List Col) -> List Both
          ++  solve(A,LB) finds a particular soln of the systems \spad{AX = B}
          ++  and a basis of the associated homogeneous systems \spad{AX = 0}
          ++  where B varies in the list of column vectors LB.

        particularSolution: (M, Col) -> PartialV
          ++ particularSolution(A,B) finds a particular solution of the linear
          ++ system \spad{AX = B}.
        hasSolution?: (M, Col) -> Boolean
          ++ hasSolution?(A,B) tests if the linear system \spad{AX = B}
          ++ has a solution.
        rank        : (M, Col) -> N
          ++ rank(A,B) computes the rank of the complete matrix \spad{(A|B)}
          ++ of the linear system \spad{AX = B}.

    Capsule ==> add
      systemMatrix      : (M, Col) -> M
      aSolution         :  M -> PartialV

      -- rank theorem
      hasSolution?(A, b) == rank A = rank systemMatrix(A, b)
      systemMatrix(m, v) == horizConcat(m, -(v::M))
      rank(A, b)         == rank systemMatrix(A, b)
      particularSolution(A, b) == aSolution rowEchelon systemMatrix(A,b)

      -- m should be in row-echelon form.
      -- last column of m is -(right-hand-side of system)
      aSolution m ==
         nvar := (ncols m - 1)::N
         rk := maxRowIndex m
         while (rk >= minRowIndex m) and every?(zero?, row(m, rk))
           repeat rk := dec rk
         rk < minRowIndex m => new(nvar, 0)
         ck := minColIndex m
         while (ck < maxColIndex m) and zero? qelt(m, rk, ck) repeat
           ck := inc ck
         ck = maxColIndex m => "failed"
         sol := new(nvar, 0)$Col
         -- find leading elements of diagonal
         v := new(nvar, minRowIndex m - 1)$PrimitiveArray(Integer)
         for i in minRowIndex m .. rk repeat
           j : Integer := 0
           while zero? qelt(m, i, j+minColIndex m) repeat j := j + 1
           v.j := i
         for j in 0..nvar-1 repeat
           if v.j >= minRowIndex m then
             qsetelt!(sol, j+minIndex sol, - qelt(m, v.j, maxColIndex m))
         sol

      solve(A:M, b:Col) ==
          -- Special case for homogeneous systems.
          every?(zero?, b) => [new(ncols A, 0), nullSpace A]
          -- General case.
          m   := rowEchelon systemMatrix(A, b)
          [aSolution m,
           nullSpace subMatrix(m, minRowIndex m, maxRowIndex m,
                                      minColIndex m, maxColIndex m - 1)]

      solve(A:M, l:List Col) ==
          null l => [[new(ncols A, 0), nullSpace A]]
          nl := (sol0 := solve(A, first l)).basis
          cons(sol0,
                 [[aSolution rowEchelon systemMatrix(A, b), nl]
                                                       for b in rest l])

@
\section{package LSMP1 LinearSystemMatrixPackage1}
<<package LSMP1 LinearSystemMatrixPackage1>>=
)abbrev package LSMP1 LinearSystemMatrixPackage1
++ Author: R. Sutor
++ Date Created: June, 1994
++ Date Last Updated:
++ Basic Functions: solve, particularSolution, hasSolution?, rank
++ Related Constructors: LinearSystemMatrixPackage
++ Also See:
++ AMS Classifications:
++ Keywords: solve
++ References:
++ Description:
++ This package solves linear system in the matrix form \spad{AX = B}.
++ It is essentially a particular instantiation of the package
++ \spadtype{LinearSystemMatrixPackage} for Matrix and Vector. This
++ package's existence makes it easier to use \spadfun{solve} in the
++ AXIOM interpreter.

LinearSystemMatrixPackage1(F): Cat == Capsule where
    F: Field
    Row      ==> Vector F
    Col      ==> Vector F
    M        ==> Matrix(F)
    LL       ==> List List F

    N        ==> NonNegativeInteger
    PartialV ==> Union(Col, "failed")
    Both     ==> Record(particular: PartialV, basis: List Col)
    LSMP     ==> LinearSystemMatrixPackage(F, Row, Col, M)

    Cat ==> with
        solve       : (M, Col) -> Both
          ++  solve(A,B) finds a particular solution of the system \spad{AX = B}
          ++  and a basis of the associated homogeneous system \spad{AX = 0}.
        solve       : (LL, Col) -> Both
          ++  solve(A,B) finds a particular solution of the system \spad{AX = B}
          ++  and a basis of the associated homogeneous system \spad{AX = 0}.
        solve       : (M, List Col) -> List Both
          ++  solve(A,LB) finds a particular soln of the systems \spad{AX = B}
          ++  and a basis of the associated homogeneous systems \spad{AX = 0}
          ++  where B varies in the list of column vectors LB.
        solve       : (LL, List Col) -> List Both
          ++  solve(A,LB) finds a particular soln of the systems \spad{AX = B}
          ++  and a basis of the associated homogeneous systems \spad{AX = 0}
          ++  where B varies in the list of column vectors LB.

        particularSolution: (M, Col) -> PartialV
          ++ particularSolution(A,B) finds a particular solution of the linear
          ++ system \spad{AX = B}.
        hasSolution?: (M, Col) -> Boolean
          ++ hasSolution?(A,B) tests if the linear system \spad{AX = B}
          ++ has a solution.
        rank        : (M, Col) -> N
          ++ rank(A,B) computes the rank of the complete matrix \spad{(A|B)}
          ++ of the linear system \spad{AX = B}.

    Capsule ==> add
        solve(m : M, c: Col): Both == solve(m,c)$LSMP
        solve(ll : LL, c: Col): Both == solve(matrix(ll)$M,c)$LSMP
        solve(m : M, l : List Col): List Both == solve(m, l)$LSMP
        solve(ll : LL, l : List Col): List Both == solve(matrix(ll)$M, l)$LSMP
        particularSolution (m : M, c : Col): PartialV == particularSolution(m, c)$LSMP
        hasSolution?(m :M, c : Col): Boolean == hasSolution?(m, c)$LSMP
        rank(m : M, c : Col): N == rank(m, c)$LSMP

@
\section{package LSPP LinearSystemPolynomialPackage}
<<package LSPP LinearSystemPolynomialPackage>>=
)abbrev package LSPP LinearSystemPolynomialPackage
++ Author:  P.Gianni
++ Date Created: Summer 1985
++ Date Last Updated: Summer 1993
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References: SystemSolvePackage
++ Description:
++ this package finds the solutions of linear systems presented as a
++ list of polynomials.

LinearSystemPolynomialPackage(R, E, OV, P): Cat == Capsule where
    R          :   IntegralDomain
    OV         :   OrderedSet
    E          :   OrderedAbelianMonoidSup
    P          :   PolynomialCategory(R,E,OV)

    F        ==> Fraction P
    NNI      ==> NonNegativeInteger
    V        ==> Vector
    M        ==> Matrix
    Soln     ==> Record(particular: Union(V F, "failed"), basis: List V F)

    Cat == with
        linSolve:  (List P, List OV) -> Soln
          ++ linSolve(lp,lvar) finds the solutions of the linear system
          ++ of polynomials lp = 0 with respect to the list of symbols lvar.

    Capsule == add

                        ---- Local Functions ----

        poly2vect:    (P,     List OV)    -> Record(coefvec: V F, reductum: F)
        intoMatrix:   (List P,   List OV) -> Record(mat: M F, vec: V F)


        poly2vect(p : P, vs : List OV) : Record(coefvec: V F, reductum: F) ==
            coefs := new(#vs, 0)$(V F)
            for v in vs for i in 1.. while p ~= 0 repeat
              u := univariate(p, v)
              degree u = 0 => "next v"
              coefs.i := (c := leadingCoefficient u)::F
              p := p - monomial(c,v, 1)
            [coefs, p :: F]

        intoMatrix(ps : List P, vs : List OV ) : Record(mat: M F, vec: V F) ==
            m := zero(#ps, #vs)$M(F)
            v := new(#ps, 0)$V(F)
            for p in ps for i in 1.. repeat
                totalDegree(p,vs) > 1 => error "The system is not linear"
                r   := poly2vect(p,vs)
                m:=setRow!(m,i,r.coefvec)
                v.i := - r.reductum
            [m, v]

        linSolve(ps, vs) ==
            r := intoMatrix(ps, vs)
            solve(r.mat, r.vec)$LinearSystemMatrixPackage(F,V F,V F,M 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 LSMP LinearSystemMatrixPackage>>
<<package LSMP1 LinearSystemMatrixPackage1>>
<<package LSPP LinearSystemPolynomialPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}