Project Euler 182

プロジェクトオイラー
http://projecteuler.net/index.php

Q182.
p = 1009、q = 3643のとき、n = pq、φ(n) = (p - 1)(q - 1)として、1 < e < φ(n)かつ(e, φ(n)) = 1のeについて、m^e ≡ m(mod n)となる0 ≤ m < nのmの個数を得たとき、最小のmの個数を持つeの総和を求めよ。

eではなく、mについて考える。(m, n) = 1のとき、m^e ≡ 1(mod n)となるeがあるから、そのような最小のeを求めればよい。そうするとeの倍数についても、この剰余式が成り立ち、eの倍数以外では成り立たない。
m^e ≡ 1(mod n)は、m^e ≡ 1(mod p)かつm^e ≡ 1(mod q)と同値。
フェルマーの小定理（この問題では間接的に使っているが、この拡張がオイラーの定理）からm^p-1 ≡ 1(mod p)より、eはp-1の約数。m^e ≡ 1(mod p)が成り立つ最小のeをe_pとすると、m^e ≡ 1(mod n)が成り立つ最小のeはe_pとe_qの最大公約数となる。
(m + kp)^e ≡ m(mod p)だから、e_p(m)は1 ≤ e_p < pについて計算しておけばよい。
eを求めた後、eの倍数についてのunconcealed messagesをインクリメントするが、そうではなく、同じeが多くあるから、同じeの数を数えておいて、まとめてunconcealed messagesを増加させれば速くなる。

import fractions

from itertools import count, ifilter
def mul(a, b, n):

    return a * b % n
def pow(m, e, n):

    if e == 1:

        return m

    

    p = pow(m, e / 2, n)

    if e & 1:

        return mul(m, mul(p, p, n), n)

    else:

        return mul(p, p, n)
def num_unconcealed_mes(e, n):

    counter = 0

    for m in xrange(1, n):

        if pow(m, e, n) == m:

            counter += 1

    return counter
def find_factor(n, p):

    while p * p <= n:

        if n % p == 0:

            return p

        if p == 2:

            p += 1

        else:

            p += 2

    return 0
def calc_divisors(n, p = 2):

    if n == 1:

        return [ 1 ]

    

    p = find_factor(n, p)

    if p:

        a = [ ]

        e = 0

        while n % p == 0:

            n /= p

            e += 1

        b = calc_divisors(n, p + 1)

        p_pow = 1

        for i in range(e + 1):

            a.extend(map(lambda x: p_pow * x, b))

            p_pow *= p

        a.sort()

        return a

    else:

        return [ 1, n ]
def LCM(m, n):

    return m * n / fractions.gcd(m, n)
def calc_cycle_core(m, p, divs):

    if m == 0:

        return 2

    

    for e in divs:

        if e == 1:

            continue

        if pow(m, e, p) == 1:

            return e + 1
def calc_cycle(m, p, q, ep, eq):

    if m % p == 0:

        return calc_cycle_core(m, q, divs_q)

    elif m % q == 0:

        return calc_cycle_core(m, p, divs_p)

    else:

        return LCM(ep[m%p], eq[m%q]) + 1
def prime_cycle(p):

    a = [ 1, 1 ]

    divs = calc_divisors(p - 1)

    for m in xrange(2, p):

        for e in divs:

            if e == 1:

                continue

            if pow(m, e, p) == 1:

                a.append(e)

                break

    return a
def calc_num(p, q):

    n = p * q

    phi = (p - 1) * (q - 1)

    ep = prime_cycle(p)

    eq = prime_cycle(q)

    a = [ 2 ] * phi

    dic_e = { }

    for m in xrange(2, n):

        e = calc_cycle(m, p, q, ep, eq)

        if e in dic_e:

            dic_e[e] += 1

        else:

            dic_e[e] = 1

    

    for e, num in dic_e.iteritems():

        for k in xrange(e, phi, e - 1):

            a[k] += num

    

    for e in xrange(2, phi):

        if fractions.gcd(e, phi) != 1:

            a[e] = n

    a[0] = a[1] = n

    return a
p = 1009

q = 3643

divs_p = calc_divisors(p - 1)

divs_q = calc_divisors(q - 1)

a = calc_num(p, q)

mn = reduce(min, a)

print sum(ifilter(lambda i: a[i] == mn, xrange(len(a))))