```-----------------------------------------------------------------------------
-- |
-- Module      :  Polynomial.Roots
--
-- 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

{-

Jan Skibinski, <a href="http://www.numeric-quest.com/news/">
Numeric Quest Inc.</a>, Huntsville, Ontario, Canada

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  =

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)

-----------------------------------------------------------------------------
--
--
--
-- Email:
--
--      jans@numeric-quest.com
--
--
--	GNU General Public License, GPL
--
-----------------------------------------------------------------------------
```