My Weblog

Blog about programming and math

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 n=
	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 
		q_1=div (n-p_1^2) q_0
		a_1=div (d+p_1) q_1
    in [a|(a,_,_)<-lst]

pellsolver n=
	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
        (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)	

	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

February 11, 2011 - Posted by | Programming

No comments yet.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your 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: