Project Euler 4（2） - inamori’s diary

Project Eulerの問題は、反対側から解いたほうがよいことが多いです。

この問題は「3桁の整数同士の積で最も大きい回文数を求めよ」ですが、これを「回文数で3桁の整数同士の積になる最も大きいものを求めよ」と読み替えます。そうすると、大きな回文数から順に調べていけばよいです。すなわち、999999,998899,997799,…と調べていき、3桁の数の積になるものが見つかれば、それが答えです。

3桁の整数同士の積になるかを判定するには、まず素因数分解をしなければなりません。そのために、あらかじめ3桁までの素数をエラトステネスのふるいで列挙し、個々に素因数分解します。まとまっていればエラトステネスのふるいを使って素因数分解できるのですが、飛び飛びの数なのでその手も使えません（10個はまとまるのだが、面倒なのでやる気がしない）。

例えば、24 = 2³ • 3なら、(2, 3, 3, 1)としておきます（本当は、( (2, 3), (3, 1) )としたいが、Pythonだと遅い）。こうしておくと簡単なコードで約数を次々と出すことが出来ます。



def gen_divisors(a, k = 0):

    if k == len(a):

        yield 1

    else:

        e = a[k+1]

        p = a[k]

        for d in gen_divisors(a, k + 2):

            for i in xrange(e + 1):

                yield d * p ** i

しかし、これでも前回より4倍ほど遅くなりました。
偶数桁の回文数は11で割り切れるので、これを使うと前回の方法はさらに速くなります。6〜8桁は次のようになります。

999000000999 = 999999 * 999001
99956644665999 = 9997647 * 9998017
9999000000009999 = 99999999 * 99990001

以下は、今回の方法のコードです。



from itertools import ifilter
def gen_digits(n):

    while n:

        yield n % 10

        n /= 10
def gen_palindrome():

    for m in xrange(10 ** N - 1, 10 ** (N - 1) - 1, -1):

        n = m

        for d in gen_digits(m):

            n = n * 10 + d

        yield n
def gen_divisors(a, k = 0):

    if k == len(a):

        yield 1

    else:

        e = a[k+1]

        p = a[k]

        for d in gen_divisors(a, k + 2):

            for i in xrange(e + 1):

                yield d * p ** i
def get_num_digits(n):

    m = 0

    while n:

        n /= 10

        m += 1

    return m
def get_exp(n, p):

    e = 0

    while n % p == 0:

        e += 1

        n /= p

    return e, n
def factorize(n):

    f = [ ]

    for p in primes:

        if p * p > n:

            break

        e, n = get_exp(n, p)

        if e:

            f.append(p)

            f.append(e)

    if n != 1:

        f.append(n)

        f.append(1)

    return f
def make_primes(max_p):

    if max_p <= 2:

        return [ 2 ]

    else:

        a = [ True ] * (max_p + 1)

        sqrt_max_p = int((max_p + 0.5) ** 0.5)

        primes = make_primes(sqrt_max_p)

        for p in primes:

            for k in xrange(p, max_p + 1, p):

                a[k] = False

        

        for k in xrange(sqrt_max_p + 1, max_p + 1):

            if a[k]:

                primes.append(k)

        return primes
def is_divisible_by_same_num_digits(n):

    f = factorize(n)

    for d1 in gen_divisors(f):

        d2 = n / d1

        if get_num_digits(d1) == N and get_num_digits(d2) == N:

            print d1, d2

            return True

    return False
N = 3

primes = make_primes(10 ** N)

for n in ifilter(is_divisible_by_same_num_digits, gen_palindrome()):

    print n

    break