[C++/Python] Project Euler #66 - Diophantine equation

This problem has a difficulty level of 25%, but without mathematical background, it’s extremely hard to solve.

If you try to just follow the formulas mechanically, this problem is very tough to crack.

In this type of indeterminate equation where the solution is not fixed, finding integer solutions is something that, apart from brute-force substitution, would be unknown to someone who has only studied middle and high school math.

The equation \(x^2 - D y^2 = 1\) is a Diophantine indeterminate equation. Depending on the value of D, there may be infinitely many solutions or not at all.

If D is a perfect square, the equation simplifies to:

\[x^2 = 1^2 + D y^2\]

But in this form, finding Pythagorean-like integer solutions that start with 1 is impossible, so the only integer solution is x = 1, y = 0.

However, if D is not a perfect square, there can be infinitely many integer solutions.

---

Diopantine's equation

In this particular problem, you are asked to find D. The answer is a value less than 1000, but finding the maximal solution may require dealing with very large numbers.

---

Here’s the source code I wrote as a reference. Note that it requires a Big Integer library.

//------------------------------------------------
//    Project Euler #66 - Diophantine equation
//        - by Aubrey Choi
//        - created at 2015-04-10
//------------------------------------------------
#include <stdio.h>
#include <memory.h>
#include <math.h>
#include "NxBigInt.h"

#define    LIMIT    1000

NxBigInt qde(int d);

int main()
{
    NxBigInt max = 0;
    int saved = 0;
    for(int d = 2 ; d <= LIMIT ; d++ )
    {
        NxBigInt x = qde(d);
        if( max < x ) { max = x; saved = d; }
    }
    printf("ans = %d\n", saved);
}

NxBigInt qde(int d)
{
    double r = sqrt(d);
    int q = 1, s = r;
    if( s*s == d ) return 0;
    NxBigInt x(r), y(1), h(1), k(0);
    while( x*x - d*y*y != 1 )
    {
        q = (d-s*s)/q;
        int a = (r+s)/q;
        s = a*q-s;

        NxBigInt ht = x, kt = y;
        x = a*x+h; h = ht;
        y = a*y+k; k = kt;
    }
    return x;
}

 

Python version:

"""
//-------------------------------------------------------------------
//    Project Euler #66 - Diophantine equation
//        - by Aubrey Choi
//        - created at 2015-04-10
//-------------------------------------------------------------------
"""
import math

def qde(d):
    r = math.sqrt(d)
    x = int(r)
    if x*x == d: return 0
    y = 1
    h = 1
    k = 0
    a = x
    q = 1
    s = x;
    while x*x - d*y*y != 1:
        q = (d-s*s)/q
        a = int((r+s)/q)
        s = a*q-s

        ht = x
        kt = y
        x = a*x+h
        h = ht
        y = a*y+k
        k = kt

    return x

max = 0
saved = 0
for d in range(1, 1001):
    x = qde(d)
    if x > max: 
        max = x
        saved = d
print 'ans = ', max, '(', saved, ')'