Project Euler 26

Problem 26
(略)1/dが10進で最も長い循環になるd < 1000を求めよ。
http://projecteuler.net/index.php?section=problems&id=26

例えば、

1/7 = 0.142857142857…

循環長は6となります。
これは割り算を進めていけばわかります。10を7で割ると余り3、30を7で割ると余り2、以下、6→4→5→1となって元に戻ります。余りが前と同じになるまでが周期です。


from itertools import count

def cycle(d):
a = [ 0 ] * d
m = 1
for n in count(1):
a[m] = 1
m = m * 10 % d
if a[m]:
return n
elif m == 0:
return 0

N = 10 ** 3
print max(map(cycle, range(2, N)))

しかし、数学の力を使うと1回ずつ割り算しなくてもよいことがわかります。
さきほどの1/7の計算手順をもう一度違った形で書くと、

10 = 7 * 1 + 3
30 = 7 * 4 + 2
20 = 7 * 2 + 6
60 = 7 * 8 + 4
40 = 7 * 5 + 5
50 = 7 * 7 + 1

これを次のように変形すると、


1000000 = 7 * 100000 + 300000
300000 = 7 * 40000 + 20000
20000 = 7 * 2000 + 6000
6000 = 7 * 800 + 400
400 = 7 * 50 + 50
50 = 7 * 7 + 1

1000000を7で割ると1余ることになります。また、

10n ≡ 1(mod 7)

となる最小の自然数nが6ということです。すなわち、1/dの循環の周期を求めるというのは、

10n ≡ 1(mod d)

を満たす最小の自然数nを求める、というのと同じことです。
また、dが2か5の因数があるときはそれを除けばよいです。例えば、28なら7にして、100/28 = 25/7 = 3 + 4/7だから、

4 • 10n ≡ 4(mod 7)

4と7は互いに素だから、

10n ≡ 1(mod 7)

となる最小の自然数nを求める問題になります。
さて、これはオイラーの定理と同じ形をしています。

aφ(n) ≡ 1 (mod n)
anは互いに素

φオイラーのφ関数です。

10φ(d) ≡ 1 (mod d)

は確かに成り立つのですが、φ(d)は

10n ≡ 1 (mod d)

を満たす最小のnとは限りません。それは、φ(d)の約数になります。これは数学的考察で簡単にわかりますが、ここでは実際に確かめてみましょう。

d d' φ(d') 周期
21 21 12 6
22 11 10 2
23 23 22 22
24 3 2 1
25 1 1 0
26 13 12 6
27 27 18 3
28 7 6 6
29 29 28 28

ここで、d'dから2と5を取り除いた自然数です。周期がφ(d')の約数になっていることが確認できます。
よって、小さい約数から順に10nを計算して、剰余が1になるnが出てきたら、それが最小のnです。べき乗の計算はバイナリ法を使えば非常に速くなります。これを大きいdから調べていって、そこまでの最長の周期がd-1になっていれば終わりです。
本当はn素数のときに循環長が最長になることは容易に想像できるので、そのときだけ計算すればよいのですが、これでも十分に速いのでここまでにしておきます。



from itertools import count, takewhile
import time

def gen_prime_candidates():
yield 2
yield 3
n = 0
while True:
yield n + 5
yield n + 7
n += 6

def calc_exp(n, p):
e = 0
while n % p == 0:
e += 1
n /= p
return e, n

# 60 = 2^2*3*5 -> (2, 3, 5)
def factorize(n):
f = ()
for p in takewhile(lambda p: p * p <= n, gen_prime_candidates()):
e, n = calc_exp(n, p)
if e:
f = f + (p, )
if n > 1:
f = f + (n, )
return f

def remove_25(n):
e, n = calc_exp(n, 2)
e, n = calc_exp(n, 5)
return n

def phi(n, f):
return reduce(lambda x, p: x * (p - 1) / p, f, n)

def gen_divisors(n, f, k = 0):
if k == len(f):
yield 1
else:
p = f[k]
m = 1
while True:
for n1 in gen_divisors(n, f, k + 1):
yield m * n1
if n % p:
break
m *= p
n /= p

# n^e % d
def pow_mod(n, e, d):
if e == 0:
return 1
m = pow(n, e >> 1, d)
p = m * m % d
if e & 1:
p = p * n % d
return p

def cycle(d):
m = remove_25(d)
if m == 1:
return 0
f = factorize(m)
ph = phi(m, f)
f2 = factorize(ph)
a = [ e for e in gen_divisors(ph, f2) ]
a.sort()
for e in a:
if pow_mod(10, e, d) == 1:
return e

t = time.clock()
N = 10 ** 7
max_cycle = 0
for d in range(N - 1, 1, -1):
max_cycle = max(max_cycle, cycle(d))
if max_cycle >= d - 1:
break
print max_cycle
print time.clock() - t