----------------------------------------------------------------------------- -- | -- Module : Polynomial.Roots -- Copyright : (c) 1998 Numeric Quest Inc., All rights reserved -- License : GPL -- -- Maintainer : m.p.donadio@ieee.org -- Stability : experimental -- Portability : portable -- -- Root finder using Laguerre's method -- ----------------------------------------------------------------------------- -- This file was sucked out of the Wayback Machine at www.archive.org. -- This was orginally a HTML files containing literate Haskell. It has -- been modified to use the Polynomial library, and Haddock style comments -- have been added. As much as the original formatting has been retained -- as possible. --mpd -- Original comments are below {- Literate Haskell module <i>Roots.lhs</i> Jan Skibinski, <a href="http://www.numeric-quest.com/news/"> Numeric Quest Inc.</a>, Huntsville, Ontario, Canada 1998.09.05, last modified 1998.09.24 This module implements <i>Laguerre's</i> method for finding complex roots of polynomials. According to [1], it <i> is by far the most straightforward of these sure-fire methods. It does require that you perform complex arithmetic (even while converging to real roots), but it is guaranteed to converge to a root from any starting point. In some instances the complex arithmetic is no disadvantage, since the polynomial itself may have complex coefficients. </i> [1] Numerical Recipes in Pascal, W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Cambridge University Press, ISBN 0-521-37516-9 See also some other variations of the same book by the same authors: Numerical Recipes in C, Fortran, etc. I just happen to own [1], although I have never programmed in Pascal. :-) Example To solve the equation x^2 - 3 x + 2 = 0 form the list of coefficients [2, -3, 1] (notice the reverse order of coefficients) and execute roots 1.0e-6 300 [2,-3, 1] -- where -- 1.0e-6 is a required accuracy -- 300 is a count of permitted iterations -- (You set it to small number just in case you -- do not trust the algorithm. But if you do, -- then set it to something big, say 300) The answer is [2.0 :+ 0.0, 1.0 :+ 0.0]; that is, both roots are real and equal to 2 and 1: x^2 - 3 x + 2 = (x - 2) (x - 1) = 0 -} module Polynomial.Roots (roots) where import Data.Complex import Polynomial.Basic -- * Functions -- | Root finder using Laguerre's method roots :: RealFloat a => a -- ^ epsilon -> Int -- ^ iteration limit -> [Complex a] -- ^ the polynomial -> [Complex a] -- ^ the roots roots eps0 count0 as0 = -- -- List of complex roots of a polynomial -- a0 + a1*x + a2*x^2... -- represented by the list as=[a0,a1,a2...] -- where -- eps is a desired accuracy -- count is a maximum count of iterations allowed -- Require: list 'as' must have at least two elements -- and the last element must not be zero roots' eps0 count0 as0 [] where roots' eps count as xs | length as <= 2 = x:xs | otherwise = roots' eps count (deflate x bs [last as]) (x:xs) where x = laguerre eps count as 0 bs = drop 1 (reverse (drop 1 as)) deflate z bs' cs | bs' == [] = cs | otherwise = deflate z (tail bs') (((head bs')+z*(head cs)):cs) laguerre :: RealFloat a => a -> Int -> [Complex a] -> Complex a -> Complex a laguerre eps0 count as0 x0 -- -- One of the roots of the polynomial 'as', -- where -- eps is a desired accuracy -- count is a maximum count of iterations allowed -- x is initial guess of the root -- This method is due to Laguerre. -- | count <= 0 = x0 | magnitude (x0 - x0') < eps0 = x0' | otherwise = laguerre eps0 (count - 1) as0 x0' where x0' = laguerre2 eps0 as0 as0' as0'' x0 as0' = polyderiv as0 as0'' = polyderiv as0' laguerre2 eps as as' as'' x -- One iteration step | magnitude b < eps = x | magnitude gp < magnitude gm = if gm == 0 then x - 1 else x - n/gm | otherwise = if gp == 0 then x - 1 else x - n/gp where gp = g + delta gm = g - delta g = d/b delta = sqrt ((n-1)*(n*h - g2)) h = g2 - f/b b = polyeval as x d = polyeval as' x f = polyeval as'' x g2 = g^(2::Int) n = fromIntegral (length as) -- Original Copyright Notice ----------------------------------------------------------------------------- -- -- Copyright: -- -- (C) 1998 Numeric Quest Inc., All rights reserved -- -- Email: -- -- jans@numeric-quest.com -- -- License: -- -- GNU General Public License, GPL -- -----------------------------------------------------------------------------