Initial population.

1 files changed, 206 insertions, 0 deletions

diff --git a/src/algebra/bezout.spad.pamphlet b/src/algebra/bezout.spad.pamphlet
new file mode 100644
index 00000000..7f6b5400
--- /dev/null
+++ b/src/algebra/bezout.spad.pamphlet
@@ -0,0 +1,206 @@
+\documentclass{article}
+\usepackage{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}