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 



February 24, 2011 Posted by | Programming | Leave a Comment

SPOJ Quadratic Equation

Spoj QUADRATE is rather parsing problem and i don’t how to parse it using more concise code.Test your parsing skills.

#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <queue>
#include <stack>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;

int main()
	{
		int t,a,b,c,ret,tmp,k_1,k_2,k_3,sign,cnt,f=0;
		for(scanf("%d\n",&t);t>0;t--)
		{
			string s;
			a=1,b=1,c=0,ret=0;	
			getline(cin,s);
			k_1=s.find("x*x",0);
			if (k_1!=0) 
			  {
				for(int i=0;i<k_1-1;i++) ret=ret*10+(s[i]-'0');
				a=ret;
			  }

			k_2=s.find("x",k_1+3);
			if(k_2==string::npos) b=0;
			else
			 {
				ret=0;
				if(s[k_1+3]=='-') sign=-1;
				else sign=1;
				for(int i=k_1+4;i<k_2 && (s[i]>='0' && s[i]<'9');i++) ret=ret*10+(s[i]-'0');
				if(ret!=0) b=ret*sign;
				else b=sign;
			 }

			k_3=s.find('=',0);
			ret=0;
			tmp=1;cnt=k_3-1;
			for(cnt=k_3-1;s[cnt]>='0' && s[cnt]<='9';cnt--)
			 {
				ret=ret+(s[cnt]-'0')*tmp;
				tmp*=10;
			 }
			if(s[cnt]=='-') c=-ret;
			else  c=ret;
			if((b*b-4*a*c)<0) cout<<"Imaginary roots.\n";
			else if ((b*b-4*a*c)==0) cout<<"Equal roots.\n";
			else cout<<"Distinct real roots.\n";
			//cout<<a<<" "<<b<<" "<<c<<endl;
			
		}
	}

February 24, 2011 Posted by | Programming | Leave a Comment

   

Follow

Get every new post delivered to your Inbox.

Join 187 other followers