My Weblog

Blog about programming and math

Naive polynomial solver in Haskell

This polynomial solver is very naive and not efficient.I wrote this solver using binary search. Solve function takes list of coefficients and determine the bound of roots and roots are polished by newton raphson method. Looking forward to implement Jenkins-Traub algorithm and Eigenvalue algorithm to locate all the roots of polynomials.

import Data.List

polydiV::(Num a)=>a->[a]->[a]
polydiV r ( x_1:x_2:[] ) = x_1 : ( x_2 + r*x_1 ) : []
polydiV r ( x_1:x_2:xs ) = x_1 : polydiV r ( ( x_2 + r*x_1 ) : xs )


polydiffreN::[Double]->[Double]
polydiffreN ps= reverse.tail.zipWith (*) [0..].reverse $ ps


horneR::Double->[Double]->Double
horneR x ps= foldl ( \seed a-> seed*x + a ) 0  ps


rootpolY::[Double]->(Double,Double)->Double
rootpolY ps (a,b) = bsearcH 0 a b where 
	bsearcH cnt lo hi 
		| cnt>=1000 =mid
		| (horneR lo ps)*(horneR mid ps)<=0 = bsearcH (cnt+1) lo mid
		| otherwise = bsearcH (cnt+1) mid hi
		where 
		   mid= (lo+hi)/2


getrangE::Double->Double->Double->[Double]->[(Double,Double)]
getrangE a b dx ps  
	| a>b = []
	|otherwise = if (horneR a ps)*(horneR (a+dx) ps)<=0 then (a,a+2*dx):getrangE (a+2*dx) b dx ps
		     else getrangE (a+dx) b dx ps

repeatrooT::[Double]->[Double]->[Double]
repeatrooT [] _ = []
repeatrooT (x:xs) ps= tmpsolvE x ps ++ repeatrooT xs ps


tmpsolvE::Double->[Double]->[Double]
tmpsolvE _ [p] = []
tmpsolvE x ps = 
   let 
	lst=polydiV x ps
   in case (abs.last $ lst) <=1.0e-6 of 
	True->x:tmpsolvE x (init lst)
        _ -> []

newTon::[Double]->[Double]->Double->(Double,Int)
newTon p p'  x_0
   | abs ( horneR x_0  p' ) <= 1.0e-9 = ( x_0 , 0 ) --p'(x_0) is zero 
   | otherwise  = until (\(_,cnt) -> cnt>=5 ) ( \( x , cnt ) -> (  x - ( ( horneR x p ) / ( horneR x p' ) ) , cnt+1 ) ) ( x_0 , 0 )



solve::[Double]->[Double]
solve  ps=
  let 
	b=foldl  max  1.0 $  map ( abs.( /(head $ ps ) ) )  ps 
	a=(-b)  
	lst=getrangE a b 0.5 ps
	tmp=map (rootpolY ps) lst
  in  if null lst then error "Sorry no real roots found" 
      else 
	 let 
	   rt = repeatrooT tmp ps
	   p=ps
	   p'= polydiffreN ps
	   proot =  map ( newTon p p' ) rt
	 in map fst  proot 



Advertisements

February 24, 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: