bancaldo™

Project Euler problem #64

All square roots are periodic when written as continued fractions
...
It can be seen that the sequence is repeating.
For conciseness, we use the notation √23 = [4;(1,3,1,8)],
to indicate that the block (1,3,1,8) repeats indefinitely.
The first ten continued fraction representations of (irrational) square roots are:

√2=[1;(2)], period=1
√3=[1;(1,2)], period=2
√5=[2;(4)], period=1
√6=[2;(2,4)], period=2
√7=[2;(1,1,1,4)], period=4
√8=[2;(1,4)], period=2
√10=[3;(6)], period=1
√11=[3;(3,6)], period=2
√12= [3;(2,6)], period=2
√13=[3;(1,1,1,1,6)], period=5

Exactly four continued fractions, for N ≤ 13, have an odd period.

How many continued fractions for N ≤ 10000 have an odd period?

Analisi:

Solo grazie a Wikipedia sono giunto alla soluzione.
Facendo riferimento a questo link, e precisamente:
“Sometimes what is desired is finding not the numerical value of a square root,
but rather its continued fraction expansion.
The following iterative algorithm can be used for this purpose
(S is any natural number that is not a perfect square):…”

In pratica il punto di partenza dell’algoritmo è:
m0 = 0
d0 = 1
a0 = la parte intera della radice quadrata del numero analizzato

>>> a0 = int(5**0.5)
>>> a0
2

a0 in sostanza è il primo termine della notazione classica delle frazioni
continue

Poi continuiamo con l’algoritmo secondo le formule:

m1 = d0 * a0 - m0
d1 = (n - m1 **2) / d0
a1 = int((a0 + m1)/ d1)

Nota: solo nel calcolo di ax, utilizziamo sempre a0, ovvero:

m2 = d1 * a1 - m1
d2 = (n - m2 **2) / d1
a2 = int((a0 + m2)/ d2)

es. con n = 114:

m0 = 0
d0 = 1
a0 = 10 (la radice quadrata di 114 = 10.677, quindi la parte intera è 10)
notazione = [10; ...]

avanti…

m1 = d0 * a0 - m0 = 1*10 - 0 = 10
d1 = (n - m1 **2) / d0 = (114 - 100) / 14 = 14
a1 = int((a0 + m1)/ d1) = int((10 + 10)/ 14 = int(20/14) = 1
notazione = [10; 1, ...]

avanti…

m2 = d1 * a1 - m1 = 14*1 - 10 = 4
d2 = (n - m2 **2) / d1 = (114 - 16) / 14 = 7
a2 = int((a0 + m2)/ d2) = int((10 + 4)/ 7 = int(17/7) = 2
notazione = [10; 1, 2, ...]

avanti…

m3 = d2 * a2 - m2 = 7*2 - 4 = 10
d3 = (n - m3 **2) / d2 = (114 - 100) / 7 = 2
a3 = int((a0 + m3)/ d3) = int((10 + 10)/ 2 = int(20/2) = 10
notazione = [10; 1, 2, 10, ...]

avanti…

m4 = d3 * a3 - m3 = 2*10 - 10 = 10
d4 = (n - m4 **2) / d3 = (114 - 100) / 2 = 7
a4 = int((a0 + m4)/ d4) = int((10 + 10)/ 7 = int(20/7) = 2
notazione = [10; 1, 2, 10, 2, ...]

avanti…

m5 = d4 * a4 - m4 = 7*2 - 10 = 4
d5 = (n - m5 **2) / d4 = (114 - 16) / 7 = 14
a5 = int((a0 + m5)/ d5) = int((10 + 4)/ 14 = int(14/14) = 1
notazione = [10; 1, 2, 10, 2, 1, ...]

avanti…

m6 = d5 * a5 - m5 = 14*1 - 4 = 10
d6 = (n - m6 **2) / d5 = (114 - 100) / 14 = 1 ** Bingo: d = 1
a6 = int((a0 + m6)/ d6) = int((10 + 10)/ 1 = int(20/1) = 20 ** Bingo: 2 * a0
notazione = [10; 1, 2, 10, 2, 1, 20, ...]

Qui, per la soluzione del nostro problema, ci fermiamo, poichè le cifre
comincerebbero a ripetersi:

m7 = d6 * a6 - m6 = 1*20 - 10 = 10 **** Uguale ad m1
d7 = (n - m7 **2) / d6 = (114 - 100) / 1 = 14 **** Uguale a d1
a7 = int((a0 + m7)/ d7) = int((10 + 10)/ 14 = int(20/14) = 1 **** Uguale a d1
notazione = [10; (1, 2, 10, 2, 1, 20)]

Quindi, ci interromperemo quando dx = 1 o quando ax = 2 * a0 (grazie Google,
anche se proseguendo in profondità e provando qualche numero, si giungerebbe
a questa soluzione anche senza il nostro immancabile amico).

Ecco il tutto messo insieme, per la risoluzione, con python, dell’odioso problema:

import time

st = time.time()
res = 0
 
for n in range(2, 10000):
    a0 = int(n ** 0.5)
    if a0 * a0 == n: continue # saltato perche' quadrato perfetto
    a, m, d = a0, 0, 1
    period = 0
    while a != 2 * a0 or d != 1: # esco dal ciclo
        m = d * a - m
        d = (n - m * m) / d
        a = int((a0 + m) / d)
        period += 1
    if period % 2 == 1: res += 1
 
print "euler 64: {}\n elapsed time: {}sec".format(res, time.time() - st)

oppure usando una funzione:

import time

st = time.time()
 
def odd_period(n):
    a0 = int(n ** 0.5)
    if a0 * a0 == n: return False # quadrato perfetto
    a, m, d = a0, 0, 1
    period = 0
    while a != 2 * a0 or d != 1:
        m = d * a - m # es. m2 = d1 * a1 - m1
        d = (n - m * m) / d # es. d2 = (n - m2*m2)/d1
        a = int((a0 + m) / d) # es. a2 = int((a0 + m2)/d2) sempre a0!
        period += 1
    if period % 2 == 1: return True

res = len([n for n in xrange(2, 10000) if odd_period(n)])
 
print "euler 64: {}\n elapsed time: {}sec".format(res, time.time() - st)

Quest’ultima è anche più veloce…
0.21sec contro 0.32sec della prima.