\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra bezout.spad}
\author{Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package BEZOUT BezoutMatrix}
<<package BEZOUT BezoutMatrix>>=
)abbrev package BEZOUT BezoutMatrix
++ Author: Clifton J. Williamson
++ Date Created: 2 August 1988
++ Date Last Updated: 3 November 1993
++ Basic Operations: bezoutMatrix, bezoutResultant, bezoutDiscriminant
++ Related Domains
++ Also See:
++ AMS Classifiactions:
++ Keywords: Bezout matrix, resultant, discriminant
++ Examples:
++ Reference: Knuth, The Art of Computer Programming, 2nd edition,
++            Vol. 2, p. 619, problem 12.
++ Description:
++   \spadtype{BezoutMatrix} contains functions for computing resultants and
++   discriminants using Bezout matrices.

BezoutMatrix(R,UP,M,Row,Col): Exports == Implementation where
  R    : Ring
  UP   : UnivariatePolynomialCategory R
  Row  : FiniteLinearAggregate R
  Col  : FiniteLinearAggregate R
  M    : MatrixCategory(R,Row,Col)
  I  ==> Integer
  lc ==> leadingCoefficient

  Exports ==> with
    sylvesterMatrix: (UP,UP) -> M
      ++ sylvesterMatrix(p,q) returns the Sylvester matrix for the two
      ++ polynomials p and q.
    bezoutMatrix: (UP,UP) -> M
      ++ bezoutMatrix(p,q) returns the Bezout matrix for the two
      ++ polynomials p and q.

    if R has commutative("*") then
      bezoutResultant: (UP,UP) -> R
       ++ bezoutResultant(p,q) computes the resultant of the two
       ++ polynomials p and q by computing the determinant of a Bezout matrix.

      bezoutDiscriminant: UP -> R
       ++ bezoutDiscriminant(p) computes the discriminant of a polynomial p
       ++ by computing the determinant of a Bezout matrix.

  Implementation ==> add

    sylvesterMatrix(p,q) ==
      n1 := degree p; n2 := degree q; n := n1 + n2
      sylmat : M := new(n,n,0)
      minR := minRowIndex sylmat; minC := minColIndex sylmat
      maxR := maxRowIndex sylmat; maxC := maxColIndex sylmat
      p0 := p
      -- fill in coefficients of 'p'
      while not zero? p0 repeat
        coef := lc p0; deg := degree p0; p0 := reductum p0
        -- put bk = coef(p,k) in sylmat(minR + i,minC + i + (n1 - k))
        for i in 0..n2 - 1 repeat
          qsetelt!(sylmat,minR + i,minC + n1 - deg + i,coef)
      q0 := q
      -- fill in coefficients of 'q'
      while not zero? q0 repeat
        coef := lc q0; deg := degree q0; q0 := reductum q0
        for i in 0..n1-1 repeat
          qsetelt!(sylmat,minR + n2 + i,minC + n2 - deg + i,coef)
      sylmat

    bezoutMatrix(p,q) ==
    -- This function computes the Bezout matrix for 'p' and 'q'.
    -- See Knuth, The Art of Computer Programming, Vol. 2, p. 619, # 12.
    -- One must have deg(p) >= deg(q), so the arguments are reversed
    -- if this is not the case.
      n1 := degree p; n2 := degree q; n := n1 + n2
      n1 < n2 => bezoutMatrix(q,p)
      m1 : I := n1 - 1; m2 : I := n2 - 1; m : I := n - 1
      -- 'sylmat' will be a matrix consisting of the first n1 columns
      -- of the standard Sylvester matrix for 'p' and 'q'
      sylmat : M := new(n,n1,0)
      minR := minRowIndex sylmat; minC := minColIndex sylmat
      maxR := maxRowIndex sylmat; maxC := maxColIndex sylmat
      p0 := p
      -- fill in coefficients of 'p'
      while not ground? p0 repeat
        coef := lc p0; deg := degree p0; p0 := reductum p0
        -- put bk = coef(p,k) in sylmat(minR + i,minC + i + (n1 - k))
        -- for i = 0...
        -- quit when i > m2 or when i + (n1 - k) > m1, whichever happens first
        for i in 0..min(m2,deg - 1) repeat
          qsetelt!(sylmat,minR + i,minC + n1 - deg + i,coef)
      q0 := q
      -- fill in coefficients of 'q'
      while not zero? q0 repeat
        coef := lc q0; deg := degree q0; q0 := reductum q0
        -- put ak = coef(q,k) in sylmat(minR + n1 + i,minC + i + (n2 - k))
        -- for i = 0...
        -- quit when i > m1 or when i + (n2 - k) > m1, whichever happens first
        -- since n2 - k >= 0, we quit when i + (n2 - k) > m1
        for i in 0..(deg + n1 - n2 - 1) repeat
          qsetelt!(sylmat,minR + n2 + i,minC + n2 - deg + i,coef)
      -- 'bezmat' will be the 'Bezout matrix' as described in Knuth
      bezmat : M := new(n1,n1,0)
      for i in 0..m2 repeat
        -- replace A_i by (b_0 A_i + ... + b_{n_2-1-i} A_{n_2 - 1}) -
        -- (a_0 B_i + ... + a_{n_2-1-i} B_{n_2-1}), as in Knuth
        bound : I := n2 - i; q0 := q
        while not zero? q0 repeat
          deg := degree q0
          if (deg < bound) then
            -- add b_deg A_{n_2 - deg} to the new A_i
            coef := lc q0
            for k in minC..maxC repeat
              c := coef * qelt(sylmat,minR + m2 - i - deg,k) +
                          qelt(bezmat,minR + m2 - i,k)
              qsetelt!(bezmat,minR + m2 - i,k,c)
          q0 := reductum q0
        p0 := p
        while not zero? p0 repeat
          deg := degree p0
          if deg < bound then
            coef := lc p0
            -- subtract a_deg B_{n_2 - deg} from the new A_i
            for k in minC..maxC repeat
              c := -coef * qelt(sylmat,minR + m - i - deg,k) +
                           qelt(bezmat,minR + m2 - i,k)
              qsetelt!(bezmat,minR + m2 - i,k,c)
          p0 := reductum p0
      for i in n2..m1 repeat for k in minC..maxC repeat
        qsetelt!(bezmat,minR + i,k,qelt(sylmat,minR + i,k))
      bezmat

    if R has commutative("*") then

      bezoutResultant(f,g) == determinant bezoutMatrix(f,g)

      if R has IntegralDomain then

        bezoutDiscriminant f ==
          degMod4 := (degree f) rem 4
          (degMod4 = 0) or (degMod4 = 1) =>
            (bezoutResultant(f,differentiate f) exquo (lc f)) :: R
          -((bezoutResultant(f,differentiate f) exquo (lc f)) :: R)

        else

          bezoutDiscriminant f ==
            lc f = 1 =>
              degMod4 := (degree f) rem 4
              (degMod4 = 0) or (degMod4 = 1) =>
                bezoutResultant(f,differentiate f)
              -bezoutResultant(f,differentiate f)
            error "bezoutDiscriminant: leading coefficient must be 1"

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