My Weblog

Blog about programming and math

Roots of polynomial equation

This is translation of Laguerre’s method given on Wikipedia. Numerical recipes suggests that use this process as approximation of root and use newton method to polish the roots. Wiki suggests “However, given the fairly limited theoretical understanding of the algorithm, many numerical analysts are hesitant to use it as such, and prefer better understood methods such as the Jenkins-Traub method, for which more solid theory has been developed”. I got one problem related to this on SPOJ [ Polynomial Equations ] but it requires lot of parsing which seems kind of hard in Haskell.

import Data.Complex
import Data.List

--complex addition
comAdd::Complex Double->Complex Double->Complex Double
comAdd u v = ( (:+) ( a + c ) ( b + d ) ) where 
	( a :+ b ) = u
	( c :+ d ) = v

--complex substraction
comSub::Complex Double->Complex Double->Complex Double
comSub u v = ( (:+) ( a - c ) ( b - d ) ) where
        ( a :+ b ) = u
        ( c :+ d ) = v

--complex multiplication
comMult::Complex Double->Complex Double->Complex Double
comMult u v =( (:+) ( a*c-b*d ) ( a*d+b*c ) ) where 
	( a :+ b ) = u
	( c :+ d ) = v

--complex division
comDiv::Complex Double->Complex Double->Complex Double
comDiv u v = ( (:+) ( ( a*c+b*d ) / k ) ( ( b*c-a*d ) / k ) ) where 
	( a :+ b ) = u
	( c :+ d ) = v
	k = c*c + d*d

--multiplication of complex number by real number
comMultn::Double->Complex Double->Complex Double
comMultn n ( a :+ b )  = ( (:+) ( n*a ) ( n*b ) ) 


--horner method to evaluate polynomial
horNer::Complex Double->[Complex Double]->Complex Double
horNer x ps = foldl ( \a b -> comAdd ( comMult a x ) b ) ( 0.0 :+ 0.0 ) ps


--divide polynomial p(x) by (x-r)
polyDiv::Complex Double->[Complex Double]->[Complex Double]
polyDiv r [x_1] = [x_1]
polyDiv r ( x_1:x_2:[] ) = x_1 : comAdd x_2 ( comMult r x_1 ):[]
polyDiv r ( x_1:x_2:xs ) = x_1 : polyDiv r ( ( comAdd x_2 ( comMult r x_1 ) ) : xs )


--differentiation of p(x)
polyDiff::[Complex Double]->[Complex Double]
polyDiff ps = reverse.tail.zipWith comMultn [0..].reverse $ ps

--sqrt of complex number
sqrT::Complex Double ->Complex Double
sqrT ( x_0 :+ y_0 )
    | x_0<0 && ( abs y_0 ) <= 1.0e-9 = ( (:+) 0.0 ( sqrt.abs $ x_0 ) )
    | otherwise = v where 
	a=sqrt ( ( x_0 + sqrt ( x_0*x_0 + y_0*y_0 ) )/2.0 )
	b= if y_0<0 then -sqrt ( ( -x_0 + sqrt ( x_0*x_0 + y_0*y_0 ) )/2.0 ) else sqrt ( ( -x_0 + sqrt ( x_0*x_0 + y_0*y_0 ) )/2.0 )
	v = ( (:+) a b )	   


--laguerre method to compute root
soLverhelp::[Complex Double]->[Complex Double]->[Complex Double]-> Complex Double-> Complex Double -> Int -> Complex Double
soLverhelp p p' p'' x_0 n cnt 
    | cnt>=100 = x_0
    | magnitude ( horNer x_0  p )  <= 1.0e-9 = x_0   --case when x_0 is itself a root 
    | otherwise = 
	let 
	   g = comDiv ( horNer x_0 p' ) ( horNer x_0 p )
	   h = comSub ( comMult g g ) ( comDiv ( horNer x_0 p'' ) ( horNer x_0 p ) )
	   t_1 = comAdd g ( sqrT ( ( comSub n 1 ) * ( comSub ( comMult n h ) ( comMult g g ) ) ) )
	   t_2 = comSub g ( sqrT ( ( comSub n 1 ) * ( comSub ( comMult n h ) ( comMult g g ) ) ) )
	   a = comDiv n ( if magnitude t_1 > magnitude t_2 then t_1 else t_2 )
	   x_0' = comSub x_0 a
	in case magnitude ( horNer x_0'  p )  <= 1.0e-9  of 
		True ->  x_0'
		_ -> soLverhelp p p' p'' x_0' n ( cnt + 1 )


--helping function
soLve::[Complex Double]->Complex Double->Complex Double
soLve ps x_0 = soLverhelp p p' p'' x_0 n cnt where 
	p = ps
	p' = polyDiff p
	p'' = polyDiff p'
	n=( genericLength ps - 1 )
	cnt = 0

--helping function
helplaGurre::[Complex Double] -> [Complex Double]
helplaGurre [x]= []
helplaGurre ps= x:helplaGurre ( init ps') where 
	x = soLve ps (1.0 :+ 0.0 )
	ps' = polyDiv x ps

--polish the roots using newton raphson method
newTon::[Complex Double]->[Complex Double]->Complex Double->(Complex Double,Int)
newTon p p'  x_0 
   | magnitude ( horNer x_0 p' ) <= 1.0e-9 = ( x_0 , 0 ) --p'(x_0) is zero 
   | otherwise	= until (\(_,cnt) -> cnt>=5 ) ( \( x , cnt ) -> ( comSub x ( comDiv ( horNer x p ) ( horNer x p' ) ) , cnt+1 ) ) ( x_0 , 0 )


laGurre::[Complex Double]->[Complex Double]
laGurre ps = 
 let 
	p = ps
	p' = polyDiff ps
	lst_1 = helplaGurre ps
	lst_2 = map ( newTon p p'  ) lst_1
 in map fst  lst_2
Advertisements

March 1, 2011 - Posted by | Programming

No comments yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: