Project Euler 70 - inamori’s diary

Problem 70
オイラーの関数φ(n)はn以下でnと互いに素な正数の個数を決めるのに用いられる。
（中略）興味深いことに、φ(87109) = 79180は、87109が79180の数字の順番を変えただけになっている。
1 < n < 10⁷でφ(n)がnの数字の順番を変えただけになっていて、n/φ(n)が最小になるものを求めよ。
http://projecteuler.net/index.php?section=problems&id=70

エラトステネスのふるい的に10⁷までのφ(n)を求め、各n/φ(n)を求めればよいでしょう（Code1）。
しかし、これはかなり遅いです。速くするにはどうしたらよいでしょう。最小を求めたいので、小さいほうから順に順列になっているかを判定して最初にこれを満たしたものが解になる、という方法でよいでしょう。これなら速くなるはずです。しかし、小さいほうから並べるのはそんなに簡単ではありません。
まず、n/φ(n)が小さくなるのは、nが素数のときでしょう。φ(p) = p - 1だからです。しかし、すぐに分かるように1を引いて順列になるような整数はありません。
次に小さいのは、p₁p₂（p_kは素数、p₁ ≤ p₂）のときでしょう。φ(p²) = p(p - 1)、φ(p₁p₂) = (p₁ - 1)(p₂ - 1)となります。3つ以上の素数から成るN未満の整数なら、p³/φ(p³) = p / (p - 1) > N^1/3 / (N^1/3 - 1)となります。これより小さい範囲なら2つの素数から成る数を調べればよいです。
これをn/φ(n)が小さい順に並べればよいのですが、ややこしいです。1000未満で考えます。p₁は31以下です。p₁p₂を(p₁, p₂)と表記することにします。p₁が同じものをまとめます。すなわち、

[(31, 31)]
[(29, 29), (29, 31)]
[(23, 23), (23, 43), (23, 41), (23, 37), (23, 31), (23, 29)]
...

実際にはジェネレータを順にリストgに入れます。そして、その中から値を一つずつリストaに入れます。最初は、gには31の系統（g31と書くことにする）を、aには(31, 31)が入ります。aは小さい順に並べてあり、aから値を一つ取り出します。取り出した系統から次の値を補充するのですが、もうg31はこれ以上値が無いので、これは捨てて次のg29をgに入れて、g29の最初の値をaに入れます。aから値を取り出したとき、次の先頭のaの値より最初の値が小さい系統をgに入れて、その最初の値をaに入れます。
こうして順番に値を出していき、p_kがN^1/3より小さくなったらそこで打ち切り、適当な方法で全部調べます（Code2）。



# Code1

from itertools import ifilter
def sieve(max_n):

    a = range(max_n + 1)

    for p in ifilter(lambda n: a[n] == n, xrange(2, max_n / 2 + 1)):

        for k in ifilter(lambda k: a[k] == k, xrange(p * 2, max_n + 1, p)):

            a[k] = p

    return a
def div_pow(n, p):

    e = 0

    while n % p == 0:

        e += 1

        n /= p

    return n, e
def calc_phi(a):

    for n in xrange(2, len(a)):

        if a[n] == n:   # prime

            a[n] -= 1

        else:

            p = a[n]

            m, e = div_pow(n, p)

            a[n] = a[p] * p ** (e - 1) * a[m]

    return a
def gen_digits(n):

    while n:

        yield n % 10

        n /= 10
def is_permutation(m, n):

    def normalize(n):

        return reduce(lambda x, y: x * 10 + y,

                    sorted(gen_digits(n), key = lambda e: -e))

    return normalize(m) == normalize(n)
def min_f(x, y):

    return x if x[0] * y[1] <= x[1] * y[0] else y
N = 10 ** 7

phi = calc_phi(sieve(N))

print time.clock() - t

g =ifilter(lambda n: is_permutation(n, phi[n]), xrange(2, N + 1))

print reduce(min_f, map(lambda n: (n, phi[n]), g), (9, 1))



# Code2

from itertools import takewhile, ifilter, imap, count

from math import sqrt

from fractions import gcd
def last(a):

    e = -1

    for e in a:

        pass

    return e
def make_primes(max_n):

    a = [ True ] * (max_n + 1)

    for p in ifilter(lambda n: a[n], xrange(2, max_n + 1)):

        for k in xrange(p * 2, max_n + 1, p):

            a[k] = False

    return filter(lambda n: a[n], xrange(2, max_n + 1))
def is_prime(n):

    return all(imap(lambda p: n % p != 0,

                takewhile(lambda p: p * p <= n, primes)))
def gen_digits(n):

    while n:

        yield n % 10

        n /= 10
def is_permutation(m, n):

    def normalize(n):

        return reduce(lambda x, y: x * 10 + y,

                    sorted(gen_digits(n), key = lambda e: -e))

    return normalize(m) == normalize(n)
def max_prime_index(limit):

    for k in xrange(len(primes)):

        if primes[k] > limit:

            return k - 1
def gen_factors(k):

    p = primes[k]

    yield (p * p, p * p - p)

    limit = N / p

    max_l = last(takewhile(lambda l: primes[l] <= limit, count(k + 1)))

    for q in imap(lambda l: primes[l], xrange(max_l, k, -1)):

        yield (p * q, (p - 1) * (q - 1))

    yield (0, p)
def gen_ratio():

    def insert_index(a, e):

        k = 0

        for k in xrange(len(a) - 1, -1, -1):

            if a[k][0][0] * e[0][1] <= a[k][0][1] * e[0][0]:

                return k + 1

        return k

    

    def insert(a, e):

        a.insert(insert_index(a, e), e)

    

    k = max_prime_index(M)

    g = [ gen_factors(k) ]

    a = [ (g[0].next(), g[0]) ]

    while k > 0:

        f, gl = a.pop(0)

        yield f

        f1 = gl.next()

        if f1[0] == 0:

            g.pop(0)

        else:

            insert(a, (f1, gl))

        

        if len(a) == 0:

            k -= 1

            g.append(gen_factors(k))

            insert(a, (g[0].next(), g[0]))

        

        # add p if p / (p - 1) < f0

        f0, g0 = a[0]

        while True:

            p = primes[k-1]

            if p * f0[1] > (p - 1) * f0[0]:

                break

            k -= 1

            g.append(gen_factors(k))

            insert(a, (g[-1].next(), g[-1]))

        if p < L:

            break
def phi(n):

    counter = 0

    for k in ifilter(lambda k: gcd(k, n) == 1, range(1, n + 1)):

        counter += 1

    return counter
N = 10 ** 7

M = int(sqrt(N))

L = int(N ** (1. / 3) + 0.5)

primes = make_primes(int(N ** (2. / 3) + 0.5))
for f in gen_ratio():

    if is_permutation(f[0], f[1]):

        print f[0]

        exit(0)
def min_f(x, y):

    return x if x[0] * y[1] < x[1] * y[0] else y
print reduce(min_f, ifilter(lambda x: is_permutation(x[0], x[1]),

                        imap(lambda n: (n, phi(n)), xrange(2, N))))