# My Weblog

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

}
}


## Project Euler problem 268

Project Euler problem 268 is related to principal of inclusion and exclusion. I was aware of inclusion exclusion but don’t how to extend it to set of k element so i did bit of brute force and got the pattern.I counted how many time a number appearing more than once in series i.e. 2310 appear 6 times. Factor of 2310 is 2*3*5*7*11.It appeared five time because of taking 4 combinations 1)2*3*5*7 (11) 2)2*3*5*11 (7) 3)2*3*7*11 (5) 4)2*5*7*11 (3) 5) 3*5*7*11 (2) and taking 6 at a time 2*3*5*7*11.Other repetitions were 12. I searched 1 4 10 since there were extra appearing in summation and OEIS. Nice problem. Happy to solve but not satisfied. Need to study combinatorics more.

#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>
#include<gmpxx.h>
typedef mpz_class _int;
using namespace std;
_int p[25]={2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97};
_int coff[25]={1, 4, 10, 20, 35, 56, 84, 120, 165, 220, 286, 364, 455, 560, 680, 816, 969, 1140, 1330, 1540, 1771, 2024, 2300, 2600, 2925};
_int Lim=1,sum=0,n;
vector<_int> C,V;

void comB(int k,_int ret,int lev,int v)
{
if(ret>Lim) return;
if(lev==v) C.push_back(ret);
else
{
for(int i=k+1;i<25;i++)
{
ret*=p[i];
comB(i,ret,lev+1,v);
ret/=p[i];
}
}

}
void print()
{
_int cnt=0,f;
for(int i=1;i+1<V.size();i++)
if((V[i-1]!=V[i] and V[i]==V[i+1]) or (V[i]==V[i+1]) or (V[i-1]==V[i] and V[i]!=V[i+1])) cnt++,f=V[i];
else { if (cnt>1) cout<<"( "<<f<<" "<<cnt<<") "; cnt=0;}

}
int main()
{
int par=1;
for(int i=0;i<16;i++) Lim*=10;
_int tmp;
for(int i=4;i<=25;i++)
{
C.clear();tmp=0;
for(int j=0;j<25;j++) comB(j,p[j],1,i);
n=C.size();
if(n==0) break;
for(int j=0;j<n;j++)
tmp+=(Lim/C[j])*par;
//for(int j=0;j<n;j++) for(_int k=C[j];k<Lim;k+=C[j]) V.push_back(k);
par*=-1;
sum+=tmp*coff[i-4];
}
//sort(V.begin(),V.end());
//print();
//cout<<endl;
cout<<sum<<endl;

}


## Baby-step Giant-step Algorithm

Baby-step giant step algorithm is used for discrete logarithm problem. I implement this code to solve spoj problem 3105 Power Modulo Inverted but i am getting TLE.For this problem i implemented a faster extended Eucleadean algorithm given on Algorithmic Number theory Vol 1 Efficient Algorithms but still out of luck :(.

import Data.List
import Data.Bits
import Data.IntMap
--a^x=b mod p find x

modpoW::Integer->Integer->Integer->Integer
modpoW a 0 n=1
modpoW a d n= fun a d where
fun a 1 = mod a n
fun a 2 = mod (a^2) n
fun a d = let
p=fun (mod (a^2) n) (shiftR d 1)
q=if (.&.) d 1 == 1 then mod (a*p) n else p
in mod q n

extendedgcD::Integer->Integer->[Integer]
extendedgcD u v=
let
m=[[1,0],[0,1]]
n=0
in helpgcD u v n m

helpgcD::Integer->Integer->Integer->[[Integer]]->[Integer]
helpgcD u v n m
| v==0 = if u<0 then [-((-1)^n)*(m!!1!!1),-((-1)^(n+1))*(m!!0!!1),-u]
else [((-1)^n)*(m!!1!!1),((-1)^(n+1))*(m!!0!!1),u]   --ux+vy=gcd(u,v)
|otherwise =
let
q= div u v
t=[[q,1],[1,0]]
m'=mult m t
u'=v
v'=u-q*v
in helpgcD u' v' (n+1) m'

--matrix multiplicatio
mult::[[Integer]]->[[Integer]]->[[Integer]]
mult m t = [[ sum $zipWith (*) x y | y<-transpose t ] | x<-m] babysteP::Integer->Integer->Integer->Integer babysteP a b p = let m=ceiling.sqrt.fromIntegral$ p
mp=Data.IntMap.fromList [(fromIntegral $modpoW a j p,fromIntegral j)|j<-[0..(m-1)]] [u,_,_]=extendedgcD a p inv=modpoW u m p y=b in recfuN y inv 0 m p mp recfuN::Integer->Integer->Integer->Integer->Integer->IntMap Key->Integer recfuN y inv cnt m p mp | cnt>=m = -1 | otherwise = case Data.IntMap.lookup (fromIntegral y) mp of Just x -> cnt*m+(fromIntegral x) _ -> recfuN (mod (y*inv) p) inv (cnt+1) m p mp main=do line<-fmap words getLine let a=read (line!!0)::Integer p=read (line!!1)::Integer b=read (line!!2)::Integer case and [a==0,b==0,p==0] of False-> do if gcd a p /=1 then putStrLn$ snd (1,"No Solution") else print $fst(babysteP a b p,"") main _ -> return ()  February 14, 2011 Posted by | Programming | Leave a Comment ## Shank Tonelli Algorithm Shanks tonelli algorithm is for finding square roots modulo prime number. It solves $x^2=n \mod p$ where n is quadratic residue (mod p ) and p is odd prime. Implementation in Haskell. import Data.List --Calculates x^2=n (mod p) where p is prime n is quadratic residue. legendrE::Integer->Integer->Integer legendrE a p = modpoW a (div (p-1) 2) p modpoW::Integer->Integer->Integer->Integer modpoW a d n= fuN a d where fuN a 1 = mod a n fuN a 2 = mod (a^2) n fuN a d = let p=fuN (mod (a^2) n) (div d 2) q=if odd d then mod (a*p) n else p in mod q n shankS::Integer->Integer->Integer shankS n p | legendrE n p /= 1 = error "non quadratic residude" | mod p 4 == 3 = n^(div (p+1) 4) |otherwise = let q = until (\x->mod x 2 /= 0) (\x->div x 2)$ pred p
s = until (\x->q*2^x == pred p) (+1) 0
z = until (\x->legendrE x p == pred p) (+1) 2
c = modpoW z q p
r = modpoW n (div (q+1) 2) p
t = modpoW n q p
m = s
in tonellI p c r t m

tonellI::Integer->Integer->Integer->Integer->Integer->Integer
tonellI p c r t m =
if mod t p == 1 then r
else let
i =repsqR (t^2) p 1 --until (\x->modpoW t (2^x) p == 1) (+1) 1
b = modpoW c (2^(m-i-1)) p
r' = mod (r*b) p
t' = mod (t*b^2) p
c' = mod (b^2) p
m' = i
in tonellI p c' r' t' m'

repsqR t p cnt| mod t p == 1 =cnt
| otherwise = repsqR (t^2) p (cnt+1)

main::IO()
main=do
print $shankS n p  February 13, 2011 Posted by | Programming | Leave a Comment ## Project Euler Problem 94 Project Euler problem 94 is similar to problem 100 but bit tricky in statement “third differs by no more than one unit”. Let if two sides have length a and third has (2b) so according to statement $a = 2*b \pm 1$ .Doing a bit analysis you will get $b = \frac { 2 + \sqrt (1+3{k^2} )} 3$ and $\frac {(-2)+\sqrt (1+3{k^2})} 3$. The only thing that bugged me for almost hour was difference by 2. The project euler forum states that (1,1,0) and (1,1,2) are not valid. import Data.List contiNuedFraction::Integer->[Integer] contiNuedFraction n= let d=truncate.sqrt.fromIntegral$ n
lst=iterate helpFun (d,0,1)
helpFun (a_0,p_0,q_0)=(a_1,p_1,q_1) where
p_1=a_0*q_0-p_0
q_1=div (n-p_1^2) q_0
a_1=div (d+p_1) q_1
in [a|(a,_,_)<-lst]

pellsolver::Integer->[(Integer,Integer)]
pellsolver n=
let
d=truncate.sqrt.fromIntegral $n lst=takeWhile (/=2*d)$ contiNuedFraction n
len=length lst
r=take (if even len then len else 2*len) $contiNuedFraction n p_0=r!!0 p_1=(r!!0)*(r!!1)+1 q_0=1 q_1=r!!1 (pfinal,_)=foldl recFun (p_1,p_0)$ drop 2 r
(qfinal,_)=foldl recFun (q_1,q_0) $drop 2 r where recFun (p_1,p_0) a=(a*p_1+p_0,p_1) in iterate (helpFun (pfinal,qfinal))$ (pfinal,qfinal) where
helpFun (x_1,y_1) (x_k,y_k)=(x_1*x_k+n*y_1*y_k,x_1*y_k+y_1*x_k)

solve =
let
lst_1=filter  (\(x,_)->mod (2+x) 3 ==0 ) $pellsolver 3 lst_2=[(6*p-2)|(x,_)<-lst_1,let p=div (x+2) 3] lst_3=filter (\(x,_)->mod (x-2) 3 ==0 )$ pellsolver 3
lst_4=[(6*p+2)|(x,_)<-lst_3,let p=div (x-2) 3]
in (sum.takeWhile (<=10^9) $lst_2) + (sum.takeWhile (<=10^9)$ lst_4)-2


## Project Euler Problem 100

Project Euler problem 100 has nothing to do with probability. At first i chose this problem only for probability sake and it turns out to Pell’s equation. Here is bit analysis. Let number of blue balls are b and red balls are r. We have
2*b*(b-1)=(b+r)*(b+r-1). Solving for b,we have $b= \frac {(1+2r)+\sqrt (1+8{r^2})} 2$ and it turns out $1+8{r^2} = {y^2}$. Find value of r and b for which $(b+r)>=10^{12}$ .I wrote Pell’s solver in Haskell for this problem and solved. 150th problem with Haskell. (pellsolver n) function generates infinite list of solution of equation $x^2-dy^2 = 1$. You may look MathWorld.

import Data.List

contiNuedFraction::Integer->[Integer]
contiNuedFraction n=
let
d=truncate.sqrt.fromIntegral $n lst=iterate helpFun (d,0,1) helpFun (a_0,p_0,q_0)=(a_1,p_1,q_1) where p_1=a_0*q_0-p_0 q_1=div (n-p_1^2) q_0 a_1=div (d+p_1) q_1 in [a|(a,_,_)<-lst] pellsolver::Integer->[(Integer,Integer)] pellsolver n= let d=truncate.sqrt.fromIntegral$ n
lst=takeWhile (/=2*d) $contiNuedFraction n len=length lst r=take (if even len then len else 2*len)$ contiNuedFraction n
p_0=r!!0
p_1=(r!!0)*(r!!1)+1
q_0=1
q_1=r!!1
(pfinal,_)=foldl recFun (p_1,p_0) $drop 2 r (qfinal,_)=foldl recFun (q_1,q_0)$ drop 2 r where
recFun (p_1,p_0) a=(a*p_1+p_0,p_1)
in iterate (helpFun (pfinal,qfinal)) $(pfinal,qfinal) where helpFun (x_1,y_1) (x_k,y_k)=(x_1*x_k+n*y_1*y_k,x_1*y_k+y_1*x_k) ans=let r=snd.head.dropWhile (\(_,r)->(r+(div (1+2*r+(truncate.sqrt.fromIntegral$ (1+8*r^2))) 2))<(10^12)) $pellsolver 8 b=div (1+2*r+(truncate.sqrt.fromIntegral$ (1+8*r^2))) 2
in b


## Project Euler Problem 153

Project euler problem 153 is one the best problem which i solved on Project Euler.I tried this problem in 2007 and could not solved.Once you spot the pattern its rather easy to solve.
1—->1
2—->1,1+i,1-i,2
3—->1,3
4—->1,1-i,1+i,2,2+2I,2-2I,4
5—->1,1+2i,1-2i,2+i,2-i,5
6—->1,1+i,1-i,2,3,3+3i,3-3i,6
7—->1,7
8—->1,1+i,1-i,2,2+2i,2-2i,4,4+4i,4-4i,8
9—->1,3,9
10—> 1,1+i,1-i,2,2+2i,2-2i,5,1+2i,1-2i,2+i,2-i,2+4i,2-4i,4+2i,4-2i,1+3i,1-3i,3+i,3-i
If you notice that 1+i and 1-i are divisors of all the numbers which are multiple of 2. 2+2i and 2-2i is divisors of all the number which are multiple of 4. 1+2i is divisor of all the numbers multiple of 5. Also out of these all the divisors only 1+i ,1-i, 1+2i,1-2i,2+i,2-i,1+3i,1-3i,3+i and 3-i are primitive solutions and rest can be derived from these solutions so i wrote brute force to test this algorithm.

#include<cstdio>
#include<iostream>
#include<map>
#include<algorithm>
using namespace std;
#define Lim 100000
typedef pair<int,int> PI;
multimap<int,int> M;
int real[Lim+1],ima[Lim+1];

int main()
{
long long sum=0;
for(int i=1;i<=Lim;i++) for(int j=i;j<=Lim;j+=i) real[j]+=i;
for(int a=1;a*a<=Lim;a++)
for(int b=1;b<=a;b++) if(__gcd(a,b)==1)M.insert(PI(a*a+b*b,(a==b)?(a+b):2*(a+b)));

for(multimap<int,int>:: iterator it=M.begin();it!=M.end();it++)
{
int i=(*it).first,val=(*it).second;
for(int j=1;i*j<=Lim;j++)
for(int k=i*j;k<=Lim;k+=i*j) ima[k]+=j*val;
}

for(int i=1;i<=Lim;i++) sum+=real[i]+ima[i];
cout<<sum<<endl;
}


But for 10^8, its not going work because this code heavily depends on memory. For this problem we don’t need any array and map. Here is source code which i derived for previous code.

#include<cstdio>
#include<iostream>
#include<map>
#include<vector>
#include<algorithm>
#define Lim 100000000
using namespace std;
typedef long long ll;
typedef pair<ll,ll> PI;

multimap<ll,ll> M;
int main()
{
ll sum=0;
for(int i=1;i<=Lim;i++)
{
sum+=(Lim/i)*i;

}

for(int a=1;a*a<Lim;a++)
for(int b=1;b<=a /*&& (a*a+b*b)<=Lim*/;b++)
if(__gcd(a,b)==1)
{
int i=(a*a+b*b),val=(a==b)?(a+b):2*(a+b);
for(int j=1;i*j<=Lim;j++)
{
sum+=(j*val*(Lim/(i*j)));
}
}

cout<<sum<<endl;
}


149th problem. One more to go to next level.

## Faster Rabin Miller Test in Haskell

I wrote this program to learn IO in Haskell and its bit faster then mine previous implementation [Technically that implementation is not fully owned by me. A lot of code was borrowed from Haskell site]. The Rabin Miller test can be written as pure function i but i rather prefer impurity :P. You can use filterM isProbable [k1..k2] to get primes in the range [k1,k2]. Remember filterM is not lazy. You can also have a look here.

import Random

helpMod::Integer->Integer->Integer->Integer
helpMod a d n= fun a d where
fun a 1 = mod a n
fun a 2 = mod (a^2) n
fun a d = let
p=fun (mod (a^2) n) (div d 2)
q=if odd d then mod (a*p) n else p
in mod q n

isProbable::Integer->IO Bool
isProbable n |n<=1 = return False
|n==2 = return True
|even n = return False
|otherwise	=
do
let d=until (\x->mod x 2/=0) (div 2) (n-1) --(n-1=2^s*d)
s=until (\x->d*2^x==pred n) (+1) 0  --counts the power
rabinMiller 0 n s d

rabinMiller::Integer->Integer->Integer->Integer->IO Bool
rabinMiller cnt n s d | cnt>=5= return True
| otherwise =
do
a<-getStdRandom $randomR (2,n-2) let x=helpMod a d n --let l="s= "++ show s ++" d= "++ show d ++" a= "++ show a ++" x= "++ show x --putStrLn l case x==1 of True -> rabinMiller (cnt+1) n s d _ -> if any ( ==pred n)$ take (fromIntegral s) \$ iterate (\e->mod (e^2) n)  x then rabinMiller (cnt+1) n s d
else return False

main = do
l<-getLine
k<-isProbable n
print k


February 5, 2011 Posted by | Programming | 1 Comment

## Pinger ???

I watched this movie long time ago. I posted it on blogspot but i am no more active there so i am posting it here since it was my first blog.

SUNDAY, JUNE 24, 2007
Today i watched Pinjar movie . Its one of the best movie i ever watched. Personally i feel that everyone in this movie had done marvelous job .Dr chandra prakash dwivedi fully justified the story without changing the original theme of Amrita pritam’s novel . We should be thankful to Amrita pritam who wrote such a touching story and draw the pain of girl who is fighting for her identity for whole life whether she is hindu or muslim. In the movie there is scene when manoj vajpayee invite a person to write the name Hamida on the hand of Paro and at that time Paro feels very bad and try to clean her hand but some how she could not succeed and the same name helps her when she was searching for Lajo and ultimately due to that name Hamida on Paro’s hand she succeed. Amazing justification, director had done his job without any partiality .During the whole movie there are lots of questions which float around me . Why a father and mother refused to accept her daughter? Only because they think that she lived with muslim and lots of questions will be posed by society if they accept her. Quite ridiculous. Who was responsible for Indo-Pak division ??? Why both can’t live in same country ?? Its only us who try to identify a person as Hindu or Muslim and give the same sense to our descendant. Why don’t we have some other parameters to identify a person? Is religion necessary for identification ? This movie makes me very sad and some times i was not able to watch some scene like when Urmila requests Manoj vajpayee to let her go home. When Lajo starts crying when she came to know about Paro and requests her to release . This movie always reminds me,you should be first human and then any thing else either Hindu or Muslim.

There may be some grammatical mistake since i wrote it long time ago and posting it as it is.