\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} <>= )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 x + an -- and q = bo x**n + b1 x**(n-1) + ... + b x + bn -- are such that q(-x**2) = p(x)p(-x), then -- bk := ak**2 + 2 * ((-1) * a*a + ... + -- (-1)**l * a*a) 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, an] aForth := cons(coefficient(p, k::NNI), aForth) aBack : L C := [] -- after k steps -- aBack = [ak, a, ..., 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 and aplus = a 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} <>= --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. @ <<*>>= <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}