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