Project Euler 26

PE 数学 プログラミング

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)
aとnは互いに素

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

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