Based on the problem description for Problem 66 of Project Euler I thought we had left the continued fractions for a while. Once I started reading up on the maths behind it and trying to solve the problem I got quite a lot wiser. But let’s start from the beginning and look at the problem description which reads

Consider quadratic Diophantine equations of the form:

x^{2}– Dy^{2}= 1

For example, when D=13, the minimal solution inxis 649^{2}– 13×180^{2}= 1.

It can be assumed that there are no solutions in positive integers when D is square.

By finding minimal solutions inxfor D = {2, 3, 5, 6, 7}, we obtain the following:

3^{2}– 2×2^{2}= 1

2^{2}– 3×1^{2}= 1

9^{2}– 5×4^{2}= 1

5^{2}– 6×2^{2}= 1

8^{2}– 7×3^{2}= 1

Hence, by considering minimal solutions inxfor D ≤ 7, the largestxis obtained when D=5.

Find the value of D ≤ 1000 in minimal solutions ofxfor which the largest value ofxis obtained.

## Analysing the problem

Knowing that this type of equation is called a Diophantine equation helps a whole lot, and so I started digging into the wikipedia article on the subject. It quickly turns out that this specific type of equation is known as Pell’s equation something where solution methods exists. There is a really good paper on that specific subject by Lenstra.

One of the solution methods is by calculating the convergents of . So if is the ith convergent of then it is the minimal solution to the equation.

So in Problem 64 we calculated the continued fraction of , once we have that calculating the convergents can be done as

Since the numerator in the continued fractions for square roots are always one. It can be initialised with

So now we have all that we need to generate solutions to the Pell’s equation, now we just need to iterate until a solution pops up.

## Implementing the solution

We have the maths and it involves quite a few variables, I have tried to make a fairly clean implementation, but I am not sure that I succeeded.

My C# solution looks like

int result = 0; BigInteger pmax = 0; for(int D = 2; D <= 1000; D++){ BigInteger limit = (int)Math.Sqrt(D); if (limit * limit == D) continue; BigInteger m = 0; BigInteger d = 1; BigInteger a = limit; BigInteger numm1 = 1; BigInteger num = a; BigInteger denm1 = 0; BigInteger den = 1; while(num*num - D*den*den != 1){ m = d * a - m; d = (D - m * m) / d; a = (limit + m) / d; BigInteger numm2 = numm1; numm1 = num; BigInteger denm2 = denm1; denm1 = den; num = a*numm1 + numm2; den = a * denm1 + denm2; } if (num > pmax) { pmax = num; result = D; } }

A lot of the code is used for keeping track of old values of the convergents. I could have made it as an array or list, but I think this will be faster.

When running the code we get

We get the largest x as a minimal solution for D=661 Solution took 14 ms

## Wrapping up

When I first saw this question it looked frightening to me, but once I started pounding on the problem it wasn’t that bad, especially since we have already encountered continued fractions before.

I know there is another way to solve it called Chakaravala’s method and Dreamshire uses that for solving the problem. It looks really neat but I can’t figure out how to programmatically do a certain step. If you know more about the method please let me know.

You can find the whole source code here. And as always I appreciate comments and questions regarding especially the problem but basically everything.

Nice, I did PE 66 only today because I am not doing problems in any order.

I just used your method for getting the period of the square root PE 64.

I don’t know how to explain my code to you, unfortunately 🙁

But here is the code, basically I only checked the square roots which have got odd period.

I forgot to point that expansion array holds the values of continued fraction expansion.

Hi Kristian

Again thanks for the great explanation.

But shouldnt the initialized values be d0 = 1 instead of n0 = 1? Have i missed something here?

Thanks

Anant

Hi Anant

Yes it should, it is correct now. Thanks for noticing.

/Kristian

First I tried to brute-force a solution, I do this usually just to see how bad it is 😉 Then I read it up and just like you I found this was the Pell’s equation and how to find minimal solutions. I must remember to always implement such problems with BigInteger, was quite a blow to me when I found out that longs just don’t do the trick…

Well, I have lost count on how many times I have failed to get the right answer due to overflow… So I follow you on that one.

The only thing you should note is that the BigInteger class is a lot slower for operations. So you only want to bring it out in case you have to represent big numbers.

hi Kristian,

would u please explain why we assumed n(-1)=1 and d(-1)=0 ?

To be honest, no I don’t remember why that is so right now.

There’s a lot of misleading, incomplete, or incorrect information floating around about the Chakravala method. The Chakravala implementation on Dreamshire is

notcorrect — the calculation of what Wikipedia calls “m” (or “p” in the Dreamshire implementation) is wrong and results in far too many iterations, although he does eventually arrive at the correct result.I did my best to cobble together all the varying information and produce what I believe is a correct implementation.

Implementation in python, running time 11.7 ms.

In the context of Project Euler – Problem 66, the following Diophantine (Pell’s) equation has been further examined.

x^2 – D*(y^2) = N

Where D = 661 and N = 1, 2, 3.

Calculations were done using ref. (3)

Case 1:

When D=661 and N=1, we get: (one solution)

x = 16421658242965910275055840472270471049

y = 638728478116949861246791167518480580

Case 2:

When D=661 and N=2, we get:

No solution

Case 3:

When D=661 and N=3, we get: (2 solutions)

x = 7313684380371167176724014194590252

y = 284469352886646877550334810886021

and

x = 3368

y = 131

___

Sources:

1) PE66, Kristian’s algorithm

http://www.mathblog.dk/files/euler/Problem66.cs

2) Big Integer Calculator

http://www.javascripter.net/math/calculators/100digitbigintcalculator.htm

3) Solving the diophantine equation x2 – Dy2=N, where D > 0 and not a perfect square

http://www.numbertheory.org/php/patzpos.html

4) Microsoft Visual C# 2010 Express

5) Microsoft Office Excel (2007)

I was wondering if anyone uses:

1) 5 variables instead of 3 ( c, f, x, y, k ) with gcd(fx,yk)=1 and gcd(xy,ck)=1. Composition((c_1,f_1,x_1,y_1,k_1),(c_2,f_2,x_2,y_2,k_2))= (1,1,c_1.f_1.x_1.c_2.f_2.x_2+D.y_1.y_2,c_1.f_1.x_1.y_2+c_2.f_2.x_2.y_2,c_1^2.f_1.k_1.c_2^2.f_2.k_2)

Reduction1(c,f,x,y,k)=(c,f,x/gcd(x,y),y/gcd(x,y),k/gcd(x,y)^2)

Reduction2(c,f,x,y,k)=(c,f.gcd(x,|k|),x/gcd(x,|k|),y,k/gcd(x,|k|)

Reduction3(c,f,x,y,k)=(c.gcd(f,|k|),f/gcd(f,|k|),x,y,k/gcd(f,|k|))

2) Previous “vectors” can be paired with the current “vector” if

(i) (v_current,v_chosen)ϵ{v_1,v2:vector;i_1,i_2,i_3:Positive Integer;p: Prime· |v_1.k|=p^i_1 and |v_2.k| = p^i_1 and 2i_3 – 1 ≥ min(i_1,i_2) and (v_1 * v_2).y ≡ 0 mod p^i_3 · (v_1, v_2), (v_2,v_1) }

(ii) (v_current,v_chosen) ϵ {v_1,v2 : vector; i_1,i_2,i_3 : Positive Integer; p: Prime set minus {2} · |v_1.k| = p^i_1 and |v_2.k| ϵ {2p^i_2,4p^i_2} and 2i_3 – 1 ≥ i_1 and (v_1 * v_2).y ≡ 0 mod p^i_3 · (v_1, v_2),(v_2,v_1) } or

(iii) (v_current,v_chosen) ϵ {v_1,v2 : vector; p: Prime set minus {2} · |v_1.k| = 2p and |v_2.k| ϵ {2p, p^2,2p^2,p^3} and (v_1 * v_2).y ≡ 0 mod p · (v_1, v_2), (v_2,v_1) }

3) If |k| ϵ {1,2,4} a solution can be determined directly. The solution equations can be obtained by composing (c,f,x,y,k) with itself a set number of times (for D=343 (7,1,3,1,2) needs 13 compositions). The most frequent solutions are

( i) 1 composition x=(2.f.x^2-k)/|k|,y=(2.x.y)/(c.|k|)

( ii) 5 compositions x=[2.f.x^2.(4.f.x^2-3k)^2-k^3]/|k|^3,y=2.x.y.(4.f.x^2-3k).(4.f.x^2-k)/(c.|k|^3)

(iii) 3 compositions x=[2.f.x^2.(2.f.x^2-k)^2-k^2]/(c.|k|^2)

( iv) 0 compositions x=(f.x)/(f.k)^1/2,y=y/

( v) 2 compositions x=f.x.(4.f.x^2-3.k)/[(f.k)^1/2.|k|],y=y.(4.f.x^2-k)/

Tell me if this is useful or not.