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
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;
}
}