open-axiom repository from github
\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}