Project Euler 66: Investigate the Diophantine equation x^2 − Dy^2 = 1.

Project Euler 66: Investigate the Diophantine equation x^2 − Dy^2 = 1.

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:

x2 – Dy2 = 1

For example, when D=13, the minimal solution in x is 6492 – 13×1802 = 1.

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

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

32 – 2×22 = 1
22 – 3×12 = 1
92 – 5×42 = 1
52 – 6×22 = 1
82 – 7×32 = 1

Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.

Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is 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.

Posted by Kristian


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.

BigInteger largestX=BigInteger.ZERO;
	    BigInteger tempX;
	    int D = 0;
	    int periodVal;

	    for(int i = 2; i<1000; i++)
	        // cout<<i<<endl;

	        periodVal = period(i);

	        if(periodVal % 2 == 0)continue;

	        tempX = solve(periodVal*2);    //solve and return P 2*n
	        if(tempX.compareTo(largestX) == 1)
	            largestX = tempX;
	            D = i;


	    System.out.println("Project Euler 66:\t" + D);


	int period(int n)

	    int a0 = (int) Math.sqrt(n);
	    int period = 0;
	    int d = 1;
	    int m = 0;
	    int a = a0;

	    totalExp =0;
	    expansion[totalExp++] = a;

	        m = d*a - m;
	        d = (n - m * m) / d;
	        a = (a0 + m) / d;
	        expansion[totalExp++] = a;
	    while(a != 2*a0);

	     for(int i = 1;i<period;i++)
	     expansion[i + period] = expansion[i];
	    return period;

	BigInteger solve(int n)
	     BigInteger numerators[] = new BigInteger[5000];
	     numerators[0] = BigInteger.valueOf(expansion[0]);
	     numerators[1] = numerators[0].multiply(BigInteger.valueOf(expansion[1]));
	     numerators[1] = numerators[1].add(BigInteger.ONE);
	     int i =2;
	     while(i < n) {
	     numerators[i] = numerators[i-1].multiply(BigInteger.valueOf(expansion[i])).add(numerators[i - 2]);
	     return numerators[n-1];


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?


Hi Anant

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


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 not correct — 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.

#! /usr/bin/env python
# -*- coding: utf-8 -*-

from math import sqrt

# See:
def chakravala(N):
    &quot;&quot;&quot; Return the minimal solution (x, y) for the Diophantine
    equation x^2 - Ny^2 == 1 &quot;&quot;&quot;
    m = m0 = int(round(sqrt(N)))
    a, b, k = m, 1, m*m - N
    if k == 0:
        # no solution if N is a perfect square
        return (None, None)

    while k != 1:
        # terminate if k in [-1, ±2], or k is ±4 and a or b is even
        if k == -1 or abs(k) == 2 or (abs(k) == 4 and not (a&amp;1 and b&amp;1)):
            # compose (a, b, k) with (a, b, k) and return solution
            return ((a*a + N*b*b)//abs(k), 2*a*b//abs(k))

        # find m such that: (a + b*m) % k == 0, abs(m^2 - N) is minimized
        diff = (m + m0) % abs(k)
        m_lo = m0 - diff
        m_hi = m_lo + abs(k)
        m = m_hi if abs(m_hi*m_hi - N) &lt; abs(m_lo*m_lo - N) else m_lo

        # compose (a, b, k) with (m, 1, m^2 - N)
        a, b, k = (m*a + N*b)//abs(k), (a + b*m)//abs(k), (m*m - N)//k
    return (a, b)

def solve(limit=1000):
    return max(xrange(2, limit+1), key=lambda D: chakravala(D)[0])

if __name__ == '__main__':
    print solve()
Jean-Marie Hachey

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
x = 3368
y = 131


1) PE66, Kristian’s algorithm
2) Big Integer Calculator
3) Solving the diophantine equation x2 – Dy2=N, where D > 0 and not a perfect square
4) Microsoft Visual C# 2010 Express
5) Microsoft Office Excel (2007)

Nicholas Colpitts

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)
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.

Carl Appellof

I figured out the continued fraction method pretty quickly after a quick trip to Wikipedia. 8ms for my code.
Thank you for mentioning the Chakravala method! It led me down a fascinating path to ancient Indian number theory, and I was finally able to code up a solution using that method. Thanks also to jeras for providing a more efficient method of calculating “m” in the iterations. My code just did a serial search for “m” and ended up being slower than the continued fraction method.

Leave a Reply