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

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

Written by on 22 October 2011

Topics: Project Euler

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.

11 Comments For This Post I'd Love to Hear Yours!

  1. Eman says:

    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];
  2. Eman says:

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

  3. anant says:

    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?


  4. Kristian says:

    Hi Anant

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


  5. Philipp says:

    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…

  6. Kristian says:

    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.

  7. Touhid_Shohan says:

    hi Kristian,

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

  8. Kristian says:

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

  9. jeras says:

    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()
  10. Jean-Marie Hachey says:

    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)

  11. Nicholas Colpitts says:

    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.

Leave a Comment Here's Your Chance to Be Heard!

You can use these html tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

You can use short tags like: [code language="language"][/code]

The code tag supports the following languages: as3, bash, coldfusion, c, cpp, csharp, css, delphi, diff,erlang, groovy, javascript, java, javafx, perl, php, plain, powershell, python, ruby, scala, sql, tex, vb, xml

You can use [latex]your formula[/latex] if you want to format something using latex

This site uses cookies.