Home > Project Euler, python > Problem Euler #66

Problem Euler #66


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)
Categorie:Project Euler, python
  1. Non c'è ancora nessun commento.
  1. No trackbacks yet.

Lascia un commento

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione / Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione / Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione / Modifica )

Google+ photo

Stai commentando usando il tuo account Google+. Chiudi sessione / Modifica )

Connessione a %s...

%d blogger cliccano Mi Piace per questo: