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

No comments yet.

Advertisements

## Leave a Reply