-----------------------------------------------------------------------------
-- |
-- 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 :: forall a. RealFloat a => a -> Int -> [Complex a] -> [Complex a]
roots a
eps0 Int
count0 [Complex a]
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
	forall {t}.
RealFloat t =>
t -> Int -> [Complex t] -> [Complex t] -> [Complex t]
roots' a
eps0 Int
count0 [Complex a]
as0 []
	where
	    roots' :: t -> Int -> [Complex t] -> [Complex t] -> [Complex t]
roots' t
eps Int
count [Complex t]
as [Complex t]
xs
	        | forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex t]
as forall a. Ord a => a -> a -> Bool
<= Int
2  = Complex t
xforall a. a -> [a] -> [a]
:[Complex t]
xs
	        | Bool
otherwise       =
                 t -> Int -> [Complex t] -> [Complex t] -> [Complex t]
roots' t
eps Int
count (forall {a}. (Eq a, Num a) => a -> [a] -> [a] -> [a]
deflate Complex t
x [Complex t]
bs [forall a. [a] -> a
last [Complex t]
as]) (Complex t
xforall a. a -> [a] -> [a]
:[Complex t]
xs)
	        where
	            x :: Complex t
x  = forall a.
RealFloat a =>
a -> Int -> [Complex a] -> Complex a -> Complex a
laguerre t
eps Int
count [Complex t]
as Complex t
0
	            bs :: [Complex t]
bs = forall a. Int -> [a] -> [a]
drop Int
1 (forall a. [a] -> [a]
reverse (forall a. Int -> [a] -> [a]
drop Int
1 [Complex t]
as))
	            deflate :: a -> [a] -> [a] -> [a]
deflate a
z [a]
bs' [a]
cs
	                | [a]
bs' forall a. Eq a => a -> a -> Bool
== []  = [a]
cs
		        | Bool
otherwise  =
                         a -> [a] -> [a] -> [a]
deflate a
z (forall a. [a] -> [a]
tail [a]
bs') (((forall a. [a] -> a
head [a]
bs')forall a. Num a => a -> a -> a
+a
zforall a. Num a => a -> a -> a
*(forall a. [a] -> a
head [a]
cs))forall a. a -> [a] -> [a]
:[a]
cs)


laguerre :: RealFloat a => a -> Int -> [Complex a] -> Complex a -> Complex a
laguerre :: forall a.
RealFloat a =>
a -> Int -> [Complex a] -> Complex a -> Complex a
laguerre a
eps0 Int
count [Complex a]
as0 Complex a
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.
	--
	| Int
count forall a. Ord a => a -> a -> Bool
<= Int
0	              = Complex a
x0
	| forall a. RealFloat a => Complex a -> a
magnitude (Complex a
x0 forall a. Num a => a -> a -> a
- Complex a
x0') forall a. Ord a => a -> a -> Bool
< a
eps0 = Complex a
x0'
	| Bool
otherwise                   = forall a.
RealFloat a =>
a -> Int -> [Complex a] -> Complex a -> Complex a
laguerre a
eps0 (Int
count forall a. Num a => a -> a -> a
- Int
1) [Complex a]
as0 Complex a
x0'
	where x0' :: Complex a
x0'    = forall {a}.
RealFloat a =>
a
-> [Complex a]
-> [Complex a]
-> [Complex a]
-> Complex a
-> Complex a
laguerre2 a
eps0 [Complex a]
as0 [Complex a]
as0' [Complex a]
as0'' Complex a
x0
	      as0' :: [Complex a]
as0'   = forall a. Num a => [a] -> [a]
polyderiv [Complex a]
as0
	      as0'' :: [Complex a]
as0''  = forall a. Num a => [a] -> [a]
polyderiv [Complex a]
as0'
	      laguerre2 :: a
-> [Complex a]
-> [Complex a]
-> [Complex a]
-> Complex a
-> Complex a
laguerre2 a
eps [Complex a]
as [Complex a]
as' [Complex a]
as'' Complex a
x
	        -- One iteration step
	        | forall a. RealFloat a => Complex a -> a
magnitude Complex a
b forall a. Ord a => a -> a -> Bool
< a
eps           = Complex a
x
	        | forall a. RealFloat a => Complex a -> a
magnitude Complex a
gp forall a. Ord a => a -> a -> Bool
< forall a. RealFloat a => Complex a -> a
magnitude Complex a
gm =
		    if Complex a
gm forall a. Eq a => a -> a -> Bool
== Complex a
0 then Complex a
x forall a. Num a => a -> a -> a
- Complex a
1 else Complex a
x forall a. Num a => a -> a -> a
- Complex a
nforall a. Fractional a => a -> a -> a
/Complex a
gm
	        | Bool
otherwise                   =
		    if Complex a
gp forall a. Eq a => a -> a -> Bool
== Complex a
0 then Complex a
x forall a. Num a => a -> a -> a
- Complex a
1 else Complex a
x forall a. Num a => a -> a -> a
- Complex a
nforall a. Fractional a => a -> a -> a
/Complex a
gp
	        where gp :: Complex a
gp    = Complex a
g forall a. Num a => a -> a -> a
+ Complex a
delta
		      gm :: Complex a
gm    = Complex a
g forall a. Num a => a -> a -> a
- Complex a
delta
		      g :: Complex a
g     = Complex a
dforall a. Fractional a => a -> a -> a
/Complex a
b
		      delta :: Complex a
delta = forall a. Floating a => a -> a
sqrt ((Complex a
nforall a. Num a => a -> a -> a
-Complex a
1)forall a. Num a => a -> a -> a
*(Complex a
nforall a. Num a => a -> a -> a
*Complex a
h forall a. Num a => a -> a -> a
- Complex a
g2))
		      h :: Complex a
h     = Complex a
g2 forall a. Num a => a -> a -> a
- Complex a
fforall a. Fractional a => a -> a -> a
/Complex a
b
		      b :: Complex a
b     = forall a. Num a => [a] -> a -> a
polyeval [Complex a]
as Complex a
x
		      d :: Complex a
d     = forall a. Num a => [a] -> a -> a
polyeval [Complex a]
as' Complex a
x
		      f :: Complex a
f     = forall a. Num a => [a] -> a -> a
polyeval [Complex a]
as'' Complex a
x
		      g2 :: Complex a
g2    = Complex a
gforall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int)
		      n :: Complex a
n     = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex a]
as)

-- Original Copyright Notice

-----------------------------------------------------------------------------
--
-- Copyright:
--
--	(C) 1998 Numeric Quest Inc., All rights reserved
--
-- Email:
--
--      jans@numeric-quest.com
--
-- License:
--
--	GNU General Public License, GPL
--
-----------------------------------------------------------------------------