axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / CRFP.spad

)abbrev package CRFP ComplexRootFindingPackage
++ Author: J. Grabmeier
++ Date Created: 31 January 1991
++ Date Last Updated: 12 April 1991
++ 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 digits)
++ 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, for example, 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) : SIG == CODE where
   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)

   SIG ==> with

     complexZeros : UP -> L C
       ++ complexZeros(p) tries to determine all complex zeros
       ++ of the polynomial p with accuracy given by the package
       ++ constant globalEps which you may change by setErrorBound.

     complexZeros : (UP, R) -> L C
       ++ complexZeros(p, eps) tries to determine all complex zeros
       ++ of the polynomial p with accuracy given by eps.

     divisorCascade : (UP,UP, Boolean) -> L FAE
       ++ divisorCascade(p,tp) assumes that degree of polynomial 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, the norm of the remainder polynomial.
       ++ If info is true, then information messages are issued.

     divisorCascade : (UP,UP) -> L FAE
       ++ divisorCascade(p,tp) assumes that degree of polynomial 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, 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 eps. An overall error bound
       ++ eps0 is determined and iterated tree-like calls
       ++ to pleskenSplit are used to get the factorization.
       ++ If info is true, then information messages are given.

     factor : (UP,R)  ->  FR UP
       ++ factor(p, eps) tries to factor p into linear factors
       ++ with error atmost eps. An overall error bound
       ++ eps0 is determined and iterated tree-like calls
       ++ to pleskenSplit are used to get the factorization.

     factor : UP  ->  FR UP
       ++ factor(p) tries to factor p into linear factors
       ++ with error atmost globalEps, the internal error bound,
       ++ which can be set by setErrorBound. An overall error bound
       ++ eps0 is determined and iterated tree-like calls
       ++ to 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 that this function depends on abs.

     pleskenSplit : (UP, R, Boolean)  ->  FR UP
       ++ pleskenSplit(poly,eps,info) determines a start polynomial start
       ++ by using "startPolynomial then it increases the exponent
       ++ n of start ** n mod poly to get an approximate factor of
       ++ 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 info is true, then information messages are issued.

     pleskenSplit : (UP, R)  ->  FR UP
       ++ pleskenSplit(poly, eps)  determines a start polynomial start
       ++ by using "startPolynomial then it increases the exponent
       ++ n of start ** n mod poly to get an approximate factor of
       ++ 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 errQuot.

     rootRadius : UP -> R
       ++ rootRadius(p) calculates the root radius of p with a
       ++ maximal error quotient of 1+globalEps, where
       ++ globalEps is the internal error bound, which can be
       ++ set by setErrorBound.

     schwerpunkt : UP ->  C
       ++ schwerpunkt(p) determines the 'Schwerpunkt' of the roots of the
       ++ polynomial p of degree n, the center of gravity, which is
       ++ coeffient of \spad{x**(n-1)} divided by
       ++ n times coefficient of \spad{x**n}.

     setErrorBound : R -> R
       ++ setErrorBound(eps) changes the internal error bound,
       -- by default being 10 ** (-20) to eps, if R is
       ++ by default being 10 ** (-3) to eps, if R is
       ++ a member in the category \spadtype{QuotientFieldCategory Integer}.
       ++ The internal globalDigits is set to
       ++ \em ceiling(1/r)**2*10 being 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, 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 factors which are otherwise 1.

   CODE ==> add

     Rep := ModMonic(C, UP)

     -- constants
     c : C
     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 ?**?: (R,Fraction Integer)->R
       -- we use this, otherwise we use 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 (fp = 1) 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 :=
           (degree fac = 1) => 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
       (degree p = 1) =>
         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)
       -- unit circle splitting!
       (r1 > (2::R)) and (rd < 1/(2::R)) => [startPoly,fp] 
       -- 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 (mD := minimumDegree lp.i) > 0 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. poly 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 (degree tp > 0)  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)
       (d = 1) => 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)
         (numberOfFactors split = 1) =>
           -- 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
           (degree newFactor = 1) =>
             result := result * nilFactor(newFactor,expOfFactor)
           listOfFactors:=cons([newFactor,expOfFactor],_
             listOfFactors)
       result

     -- implementation of local functions

     absC c == nthRoot(norm(c)$C,2)

     absR r ==
       r < 0 => -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 coerce: % -> 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

)fin
     -- 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.