open-axiom repository from github
\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra crfp.spad}
\author{Johannes Grabmeier}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package CRFP ComplexRootFindingPackage}
<<package CRFP ComplexRootFindingPackage>>=
)abbrev package CRFP ComplexRootFindingPackage
++ Author: J. Grabmeier
++ Date Created: 31 January 1991
++ Date Last Updated: 12 April 1991
++ Basic Operations: factor, pleskenSplit
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords: complex zeros, roots
++ References: J. Grabmeier: On Plesken's root finding algorithm,
++ in preparation
++ A. Schoenhage: The fundamental theorem of algebra in terms of computational
++ complexity, preliminary report, Univ. Tuebingen, 1982
++ Description:
++ \spadtype{ComplexRootFindingPackage} provides functions to
++ find all roots of a polynomial p over the complex number by
++ using Plesken's idea to calculate in the polynomial ring
++ modulo f and employing the Chinese Remainder Theorem.
++ In this first version, the precision (see \spadfunFrom{digits}{Float})
++ is not increased when this is necessary to
++ avoid rounding errors. Hence it is the user's responsibility to
++ increase the precision if necessary.
++ Note also, if this package is called with e.g. \spadtype{Fraction Integer},
++ the precise calculations could require a lot of time.
++ Also note that evaluating the zeros is not necessarily a good check
++ whether the result is correct: already evaluation can cause
++ rounding errors.
ComplexRootFindingPackage(R, UP): public == private where
-- R : Join(Field, OrderedRing, CharacteristicZero)
-- Float not in CharacteristicZero !|
R : Join(Field, OrderedRing)
UP : UnivariatePolynomialCategory Complex R
C ==> Complex R
FR ==> Factored
I ==> Integer
L ==> List
FAE ==> Record(factors : L UP, error : R)
NNI ==> NonNegativeInteger
OF ==> OutputForm
ICF ==> IntegerCombinatoricFunctions(I)
public ==> with
complexZeros : UP -> L C
++ complexZeros(p) tries to determine all complex zeros
++ of the polynomial p with accuracy given by the package
++ constant {\em globalEps} which you may change by
++ {\em setErrorBound}.
complexZeros : (UP, R) -> L C
++ complexZeros(p, eps) tries to determine all complex zeros
++ of the polynomial p with accuracy given by {\em eps}.
divisorCascade : (UP,UP, Boolean) -> L FAE
++ divisorCascade(p,tp) assumes that degree of polynomial {\em tp}
++ is smaller than degree of polynomial p, both monic.
++ A sequence of divisions are calculated
++ using the remainder, made monic, as divisor
++ for the the next division. The result contains also the error of the
++ factorizations, i.e. the norm of the remainder polynomial.
++ If {\em info} is {\em true}, then information messages are issued.
divisorCascade : (UP,UP) -> L FAE
++ divisorCascade(p,tp) assumes that degree of polynomial {\em tp}
++ is smaller than degree of polynomial p, both monic.
++ A sequence of divisions is calculated
++ using the remainder, made monic, as divisor
++ for the the next division. The result contains also the error of the
++ factorizations, i.e. the norm of the remainder polynomial.
factor: (UP,R,Boolean) -> FR UP
++ factor(p, eps, info) tries to factor p into linear factors
++ with error atmost {\em eps}. An overall error bound
++ {\em eps0} is determined and iterated tree-like calls
++ to {\em pleskenSplit} are used to get the factorization.
++ If {\em info} is {\em true}, then information messages are given.
factor: (UP,R) -> FR UP
++ factor(p, eps) tries to factor p into linear factors
++ with error atmost {\em eps}. An overall error bound
++ {\em eps0} is determined and iterated tree-like calls
++ to {\em pleskenSplit} are used to get the factorization.
factor: UP -> FR UP
++ factor(p) tries to factor p into linear factors
++ with error atmost {\em globalEps}, the internal error bound,
++ which can be set by {\em setErrorBound}. An overall error bound
++ {\em eps0} is determined and iterated tree-like calls
++ to {\em pleskenSplit} are used to get the factorization.
graeffe : UP -> UP
++ graeffe p determines q such that \spad{q(-z**2) = p(z)*p(-z)}.
++ Note that the roots of q are the squares of the roots of p.
norm : UP -> R
++ norm(p) determines sum of absolute values of coefficients
++ Note: this function depends on \spadfunFrom{abs}{Complex}.
pleskenSplit: (UP, R, Boolean) -> FR UP
++ pleskenSplit(poly,eps,info) determines a start polynomial {\em start}
++ by using "startPolynomial then it increases the exponent
++ n of {\em start ** n mod poly} to get an approximate factor of
++ {\em poly}, in general of degree "degree poly -1". Then a divisor
++ cascade is calculated and the best splitting is chosen, as soon
++ as the error is small enough.
--++ In a later version we plan
--++ to use the whole information to get a split into more than 2
--++ factors.
++ If {\em info} is {\em true}, then information messages are issued.
pleskenSplit: (UP, R) -> FR UP
++ pleskenSplit(poly, eps) determines a start polynomial {\em start}\
++ by using "startPolynomial then it increases the exponent
++ n of {\em start ** n mod poly} to get an approximate factor of
++ {\em poly}, in general of degree "degree poly -1". Then a divisor
++ cascade is calculated and the best splitting is chosen, as soon
++ as the error is small enough.
--++ In a later version we plan
--++ to use the whole information to get a split into more than 2
--++ factors.
reciprocalPolynomial: UP -> UP
++ reciprocalPolynomial(p) calulates a polynomial which has exactly
++ the inverses of the non-zero roots of p as roots, and the same
++ number of 0-roots.
rootRadius: (UP,R) -> R
++ rootRadius(p,errQuot) calculates the root radius of p with a
++ maximal error quotient of {\em errQuot}.
rootRadius: UP -> R
++ rootRadius(p) calculates the root radius of p with a
++ maximal error quotient of {\em 1+globalEps}, where
++ {\em globalEps} is the internal error bound, which can be
++ set by {\em setErrorBound}.
schwerpunkt: UP -> C
++ schwerpunkt(p) determines the 'Schwerpunkt' of the roots of the
++ polynomial p of degree n, i.e. the center of gravity, which is
++ {\em coeffient of \spad{x**(n-1)}} divided by
++ {\em n times coefficient of \spad{x**n}}.
setErrorBound : R -> R
++ setErrorBound(eps) changes the internal error bound,
-- by default being {\em 10 ** (-20)} to eps, if R is
++ by default being {\em 10 ** (-3)} to eps, if R is
++ a member in the category \spadtype{QuotientFieldCategory Integer}.
++ The internal {\em globalDigits} is set to
-- {\em ceiling(1/r)**2*10} being {\em 10**41} by default.
++ {\em ceiling(1/r)**2*10} being {\em 10**7} by default.
startPolynomial: UP -> Record(start: UP, factors: FR UP)
++ startPolynomial(p) uses the ideas of Schoenhage's
++ variant of Graeffe's method to construct circles which separate
++ roots to get a good start polynomial, i.e. one whose
++ image under the Chinese Remainder Isomorphism has both entries
++ of norm smaller and greater or equal to 1. In case the
++ roots are found during internal calculations.
++ The corresponding factors
++ are in {\em factors} which are otherwise 1.
private ==> add
Rep := ModMonic(C, UP)
-- constants
r : R
--globalDigits : I := 10 ** 41
globalDigits : I := 10 ** 7
globalEps : R :=
--a : R := (1000000000000000000000 :: I) :: R
a : R := (1000 :: I) :: R
1/a
emptyLine : OF := " "
dashes : OF := center "---------------------------------------------------"
dots : OF := center "..................................................."
one : R := 1$R
two : R := 2 * one
ten : R := 10 * one
eleven : R := 11 * one
weakEps := eleven/ten
--invLog2 : R := 1/log10 (2*one)
-- signatures of local functions
absC : C -> R
--
absR : R -> R
--
calculateScale : UP -> R
--
makeMonic : UP -> UP
-- 'makeMonic p' divides 'p' by the leading coefficient,
-- to guarantee new leading coefficient to be 1$R we cannot
-- simply divide the leading monomial by the leading coefficient
-- because of possible rounding errors
min: (FAE, FAE) -> FAE
-- takes factorization with smaller error
nthRoot : (R, NNI) -> R
-- nthRoot(r,n) determines an approximation to the n-th
-- root of r, if \spadtype{R} has {\em ?**?: (R,Fraction Integer)->R}
-- we use this, otherwise we use {\em approxNthRoot} via
-- \spadtype{Integer}
shift: (UP,C) -> UP
-- shift(p,c) changes p(x) into p(x+c), thereby modifying the
-- roots u_j of p to the roots (u_j - c) of shift(p,c)
scale: (UP,C) -> UP
-- scale(p,c) changes p(x) into p(cx), thereby modifying the
-- roots u_j of p to the roots ((1/c) u_j) of scale(p,c)
-- implementation of exported functions
complexZeros(p,eps) ==
--r1 : R := rootRadius(p,weakEps)
--eps0 : R = r1 * nthRoot(eps, degree p)
-- right now we are content with
eps0 : R := eps/(ten ** degree p)
facs : FR UP := factor(p,eps0)
[-coefficient(linfac.factor,0) for linfac in factors facs]
complexZeros p == complexZeros(p,globalEps)
setErrorBound r ==
r <= 0 => error "setErrorBound: need error bound greater 0"
globalEps := r
if R has QuotientFieldCategory Integer then
rd : Integer := ceiling(1/r)
globalDigits := rd * rd * 10
lof : List OF := _
["setErrorBound: internal digits set to",globalDigits::OF]
print hconcat lof
messagePrint "setErrorBound: internal error bound set to"
globalEps
pleskenSplit(poly,eps,info) ==
p := makeMonic poly
fp : FR UP
if not zero? (md := minimumDegree p) then
fp : FR UP := irreducibleFactor(monomial(1,1)$UP,md)$(FR UP)
p := p quo monomial(1,md)$UP
sP : Record(start: UP, factors: FR UP) := startPolynomial p
fp : FR UP := sP.factors
if not one? fp then
qr: Record(quotient: UP, remainder: UP):= divide(p,makeMonic expand fp)
p := qr.quotient
st := sP.start
zero? degree st => fp
-- we calculate in ModMonic(C, UP),
-- next line defines the polynomial, which is used for reducing
setPoly p
nm : R := eps
split : FAE
sR : Rep := st :: Rep
psR : Rep := sR ** (degree poly)
notFoundSplit : Boolean := true
while notFoundSplit repeat
-- if info then
-- lof : L OF := ["not successfull, new exponent:", nn::OF]
-- print hconcat lof
psR := psR * psR * sR -- exponent (2*d +1)
-- be careful, too large exponent results in rounding errors
-- tp is the first approximation of a divisor of poly:
tp : UP := lift psR
zero? degree tp =>
if info then print "we leave as we got constant factor"
nilFactor(poly,1)$(FR UP)
-- this was the case where we don't find a non-trivial factorization
-- we refine tp by repeated polynomial division and hope that
-- the norm of the remainder gets small from time to time
splits : L FAE := divisorCascade(p, makeMonic tp, info)
split := reduce(min,splits)
notFoundSplit := (eps <= split.error)
for fac in split.factors repeat
fp :=
one? degree fac => fp * nilFactor(fac,1)$(FR UP)
fp * irreducibleFactor(fac,1)$(FR UP)
fp
startPolynomial p == -- assume minimumDegree is 0
--print (p :: OF)
fp : FR UP := 1
one? degree p =>
p := makeMonic p
[p,irreducibleFactor(p,1)]
startPoly : UP := monomial(1,1)$UP
eps : R := weakEps -- 10 per cent errors allowed
r1 : R := rootRadius(p, eps)
rd : R := 1/rootRadius(reciprocalPolynomial p, eps)
(r1 > (2::R)) and (rd < 1/(2::R)) => [startPoly,fp] -- unit circle splitting!
-- otherwise the norms of the roots are too closed so we
-- take the center of gravity as new origin:
u : C := schwerpunkt p
startPoly := startPoly-monomial(u,0)
p := shift(p,-u)
-- determine new rootRadius:
r1 : R := rootRadius(p, eps)
startPoly := startPoly/(r1::C)
-- use one of the 4 points r1*zeta, where zeta is a 4th root of unity
-- as new origin, this could be changed to an arbitrary list
-- of elements of norm 1.
listOfCenters : L C := [complex(r1,0), complex(0,r1), _
complex(-r1,0), complex(0,-r1)]
lp : L UP := [shift(p,v) for v in listOfCenters]
-- next we check if one of these centers is a root
centerIsRoot : Boolean := false
for i in 1..maxIndex lp repeat
if positive? (mD := minimumDegree lp.i) then
pp : UP := monomial(1,1)-monomial(listOfCenters.i-u,0)
centerIsRoot := true
fp := fp * irreducibleFactor(pp,mD)
centerIsRoot =>
p := shift(p,u) quo expand fp
--print (p::OF)
zero? degree p => [p,fp]
sP:= startPolynomial(p)
[sP.start,fp]
-- choose the best one w.r.t. maximal quotient of norm of largest
-- root and norm of smallest root
lpr1 : L R := [rootRadius(q,eps) for q in lp]
lprd : L R := [1/rootRadius(reciprocalPolynomial q,eps) for q in lp]
-- later we should check here of an rd is smaller than globalEps
lq : L R := []
for i in 1..maxIndex lpr1 repeat
lq := cons(lpr1.i/lprd.i, lq)
--lq : L R := [(l/s)::R for l in lpr1 for s in lprd])
lq := reverse lq
po := position(reduce(max,lq),lq)
--p := lp.po
--lrr : L R := [rootRadius(p,i,1+eps) for i in 2..(degree(p)-1)]
--lrr := concat(concat(lpr1.po,lrr),lprd.po)
--lu : L R := [(lrr.i + lrr.(i+1))/2 for i in 1..(maxIndex(lrr)-1)]
[startPoly - monomial(listOfCenters.po,0),fp]
norm p ==
-- reduce(_+$R,map(absC,coefficients p))
nm : R := 0
for c in coefficients p repeat
nm := nm + absC c
nm
pleskenSplit(poly,eps) == pleskenSplit(poly,eps,false)
graeffe p ==
-- If p = ao x**n + a1 x**(n-1) + ... + a<n-1> x + an
-- and q = bo x**n + b1 x**(n-1) + ... + b<n-1> x + bn
-- are such that q(-x**2) = p(x)p(-x), then
-- bk := ak**2 + 2 * ((-1) * a<k-1>*a<k+1> + ... +
-- (-1)**l * a<l>*a<l>) where l = min(k, n-k).
-- graeffe(p) constructs q using these identities.
n : NNI := degree p
aForth : L C := []
for k in 0..n repeat -- aForth = [a0, a1, ..., a<n-1>, an]
aForth := cons(coefficient(p, k::NNI), aForth)
aBack : L C := [] -- after k steps
-- aBack = [ak, a<k-1>, ..., a1, a0]
gp : UP := 0$UP
for k in 0..n repeat
ak : C := first aForth
aForth := rest aForth
aForthCopy : L C := aForth -- we iterate over aForth and
aBackCopy : L C := aBack -- aBack but do not want to
-- destroy them
sum : C := 0
const : I := -1 -- after i steps const = (-1)**i
for aminus in aBack for aplus in aForth repeat
-- after i steps aminus = a<k-i> and aplus = a<k+i>
sum := sum + const * aminus * aplus
aForthCopy := rest aForthCopy
aBackCopy := rest aBackCopy
const := -const
gp := gp + monomial(ak*ak + 2 * sum, (n-k)::NNI)
aBack := cons(ak, aBack)
gp
rootRadius(p,errorQuotient) ==
errorQuotient <= 1$R =>
error "rootRadius: second Parameter must be greater than 1"
pp : UP := p
rho : R := calculateScale makeMonic pp
rR : R := rho
pp := makeMonic scale(pp,complex(rho,0$R))
expo : NNI := 1
d : NNI := degree p
currentError: R := nthRoot(2::R, 2)
currentError := d*20*currentError
while nthRoot(currentError, expo) >= errorQuotient repeat
-- if info then print (expo :: OF)
pp := graeffe pp
rho := calculateScale pp
expo := 2 * expo
rR := nthRoot(rho, expo) * rR
pp := makeMonic scale(pp,complex(rho,0$R))
rR
rootRadius(p) == rootRadius(p, 1+globalEps)
schwerpunkt p ==
zero? p => 0$C
zero? (d := degree p) => error _
"schwerpunkt: non-zero const. polynomial has no roots and no schwerpunkt"
-- coeffient of x**d and x**(d-1)
lC : C := coefficient(p,d) -- ~= 0
nC : C := coefficient(p,(d-1) pretend NNI)
(denom := recip ((d::I::C)*lC)) case "failed" => error "schwerpunkt: _
degree * leadingCoefficient not invertible in ring of coefficients"
- (nC*(denom::C))
reciprocalPolynomial p ==
zero? p => 0
d : NNI := degree p
md : NNI := d+minimumDegree p
lm : L UP := [monomial(coefficient(p,i),(md-i) :: NNI) for i in 0..d]
sol := reduce(_+, lm)
divisorCascade(p, tp, info) ==
lfae : L FAE := nil()
for i in 1..degree tp while positive? degree tp repeat
-- USE monicDivide !!!
qr : Record(quotient: UP, remainder: UP) := divide(p,tp)
factor1 : UP := tp
factor2 : UP := makeMonic qr.quotient
-- refinement of tp:
tp := qr.remainder
nm : R := norm tp
listOfFactors : L UP := cons(factor2,nil()$(L UP))
listOfFactors := cons(factor1,listOfFactors)
lfae := cons( [listOfFactors,nm], lfae)
if info then
--lof : L OF := [i :: OF,"-th division:"::OF]
--print center box hconcat lof
print emptyLine
lof : L OF := ["error polynomial has degree " ::OF,_
(degree tp)::OF, " and norm " :: OF, nm :: OF]
print center hconcat lof
lof : L OF := ["degrees of factors:" ::OF,_
(degree factor1)::OF," ", (degree factor2)::OF]
print center hconcat lof
if info then print emptyLine
reverse lfae
divisorCascade(p, tp) == divisorCascade(p, tp, false)
factor(poly,eps) == factor(poly,eps,false)
factor(p) == factor(p, globalEps)
factor(poly,eps,info) ==
result : FR UP := coerce monomial(leadingCoefficient poly,0)
d : NNI := degree poly
--should be
--den : R := (d::I)::R * two**(d::Integer) * norm poly
--eps0 : R := eps / den
-- for now only
eps0 : R := eps / (ten*ten)
one? d => irreducibleFactor(poly,1)$(FR UP)
listOfFactors : L Record(factor: UP,exponent: I) :=_
list [makeMonic poly,1]
if info then
lof : L OF := [dashes,dots,"list of Factors:",dots,listOfFactors::OF, _
dashes, "list of Linear Factors:", dots, result::OF, _
dots,dashes]
print vconcat lof
while not null listOfFactors repeat
p : UP := (first listOfFactors).factor
exponentOfp : I := (first listOfFactors).exponent
listOfFactors := rest listOfFactors
if info then
lof : L OF := ["just now we try to split the polynomial:",p::OF]
print vconcat lof
split : FR UP := pleskenSplit(p, eps0, info)
one? numberOfFactors split =>
-- in a later version we will change error bound and
-- accuracy here to deal this case as well
lof : L OF := ["factor: couldn't split factor",_
center(p :: OF), "with required error bound"]
print vconcat lof
result := result * nilFactor(p, exponentOfp)
-- now we got 2 good factors of p, we drop p and continue
-- with the factors, if they are not linear, or put a
-- linear factor to the result
for rec in factors(split)$(FR UP) repeat
newFactor : UP := rec.factor
expOfFactor := exponentOfp * rec.exponent
one? degree newFactor =>
result := result * nilFactor(newFactor,expOfFactor)
listOfFactors:=cons([newFactor,expOfFactor],_
listOfFactors)
result
-- implementation of local functions
absC c == nthRoot(norm(c)$C,2)
absR r ==
negative? r => -r
r
min(fae1,fae2) ==
fae2.error < fae1.error => fae2
fae1
calculateScale p ==
d := degree p
maxi :R := 0
for j in 1..d for cof in rest coefficients p repeat
-- here we need abs: R -> R
rc : R := absR real cof
ic : R := absR imag cof
locmax: R := max(rc,ic)
maxi := max( nthRoot( locmax/(binomial(d,j)$ICF::R), j), maxi)
-- Maybe I should use some type of logarithm for the following:
maxi = 0$R => error("Internal Error: scale cannot be 0")
rho :R := one
rho < maxi =>
while rho < maxi repeat rho := ten * rho
rho / ten
while maxi < rho repeat rho := rho / ten
rho = 0 => one
rho
makeMonic p ==
p = 0 => p
monomial(1,degree p)$UP + (reductum p)/(leadingCoefficient p)
scale(p, c) ==
-- eval(p,cx) is missing !!
eq : Equation UP := equation(monomial(1,1), monomial(c,1))
eval(p,eq)
-- improvement?: direct calculation of the new coefficients
shift(p,c) ==
rhs : UP := monomial(1,1) + monomial(c,0)
eq : Equation UP := equation(monomial(1,1), rhs)
eval(p,eq)
-- improvement?: direct calculation of the new coefficients
nthRoot(r,n) ==
R has RealNumberSystem => r ** (1/n)
R has QuotientFieldCategory Integer =>
den : I := approxNthRoot(globalDigits * denom r ,n)$IntegerRoots(I)
num : I := approxNthRoot(globalDigits * numer r ,n)$IntegerRoots(I)
num/den
-- the following doesn't compile
--R has CoercibleTo Fraction Integer =>
-- q : Fraction Integer := coerce(r)@Fraction(Integer)
-- den : I := approxNthRoot(globalDigits * denom q ,n)$IntegerRoots(I)
-- num : I := approxNthRoot(globalDigits * numer q ,n)$IntegerRoots(I)
-- num/den
r -- this is nonsense, perhaps a Newton iteration for x**n-r here
)if false
-- for late use:
graeffe2 p ==
-- substitute x by -x :
eq : Equation UP := equation(monomial(1,1), monomial(-1$C,1))
pp : UP := p*eval(p,eq)
gp : UP := 0$UP
while pp ~= 0 repeat
i:NNI := (degree pp) quo (2::NNI)
coef:C:=
even? i => leadingCoefficient pp
- leadingCoefficient pp
gp := gp + monomial(coef,i)
pp := reductum pp
gp
shift2(p,c) ==
d := degree p
cc : C := 1
coef := List C := [cc := c * cc for i in 1..d]
coef := cons(1,coef)
coef := [coefficient(p,i)*coef.(1+i) for i in 0..d]
res : UP := 0
for j in 0..d repeat
cc := 0
for i in j..d repeat
cc := cc + coef.i * (binomial(i,j)$ICF :: R)
res := res + monomial(cc,j)$UP
res
scale2(p,c) ==
d := degree p
cc : C := 1
coef := List C := [cc := c * cc for i in 1..d]
coef := cons(1,coef)
coef := [coefficient(p,i)*coef.(i+1) for i in 0..d]
res : UP := 0
for i in 0..d repeat res := res + monomial(coef.(i+1),i)$UP
res
scale2: (UP,C) -> UP
shift2: (UP,C) -> UP
graeffe2 : UP -> UP
++ graeffe2 p determines q such that \spad{q(-z**2) = p(z)*p(-z)}.
++ Note that the roots of q are the squares of the roots of p.
)endif
@
\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 CRFP ComplexRootFindingPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}