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

\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}