bancaldo™

Consider quadratic Diophantine equations of the form:

x^2 – Dy^2 = 1

For example, when D=13, the minimal solution in x is 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 in x for 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 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.

Analisi:

In pratica la soluzione del problema, sta nell’unione dei problemi 64 e 65.
Questo perchè da wikipedia risulta che: la soluzione all’equazione di Pell

x^2 - D*y^2 = 1

sta nella ricerca delle convergenti della frazione continua della radice
quadrata di D. Grazie alle dovute semplificazioni, risulta che i valori di x
sono niente altro che i numeratori delle convergenti, mentre i valori di y,
i denominatori.

Quindi: prima si ottiene la notazione della frazione continua della radice quadrata di D, poi
si calcolano le convergenti di tale frazione continua.
Dopodichè trovare la soluzione ovvero la prima coppia di x,y per la quale si ha
x^2 – D*y^2 = 1
è un gioco abbastanza semplice.

Unica modifica sta nel fatto che non è detto che x,y si riescano a calcolare dal
periodo ciclico segnato in notazione.
Es: per D = 13, notazione [3; 1, 1, 1, 1, 6], calcolando i convergenti per
[3; 1, 1, 1, 1, 6], non si risolve l’equazione di Pell.
Si deve quindi ripetere il periodo, diventando:
[3; 1, 1, 1, 1, 6, 1, 1, 1, 1, 6]
e così via, fino alla risoluzione dell’equazione di Pell.
Sempre per D=13 si otterranno così x=649 e y=180.

python:

import time

st = time.time()

def cont_fract(n):
    '''get the continued fraction'''
    notation = []
    a0 = int(n ** 0.5)
    notation.append(a0)
    if a0 * a0 == n:
        return notation
    a, m, d = a0, 0, 1
    while a != 2 * a0 or d != 1:
        m = d * a - m
        d = (n - m * m) / d
        a = int((a0 + m) / d)
        notation.append(a)
    return notation

def get_x_pell(notation, d):
    '''get numerator of the convergents for
       wich Pell's equation is solved'''
    hs = [0,1]
    ks = [1,0]
    i = 0
    period = notation[1:]
    for n in notation:
        h = n*hs[i+1] + hs[i]
        k = n*ks[i+1] + ks[i]
        hs.append(h)
        ks.append(k)
        i += 1
        if h**2 - d*k**2 == 1:
            return h
        else: # it multiplies the period in notation
            notation.extend(period)

xs = {}
for d in xrange(2, 1001):
    if int(d**0.5) * int(d**0.5) == d: continue
    cf = cont_fract(d) # continued fraction
    x = get_x_pell(cf, d) # get convergent numerator
    xs[d] = x

max_x = max(xs.values())
res = [d for d in xs.keys() if xs[d] == max_x][0]

print 'euler 66: {}\nelapsed time: {}sec'.format(res, time.time() - st)