Project Euler 104

http://projecteuler.net/index.php?section=problems&id=104


pandigitalかどうかの判定は、各桁の数字を用意して、それを整数の各ビットに見立てて1〜9で和を取ると、

21 + 22 + ... + 29 = 1022

となるので、これで判定できます。
下9桁の数字を出すのは、まず10^9の剰余を出して、

def digits(n):
    return [] if n == 0 else digits(n / 10) + [ n % 10 ]

とします。ただ、やはりこれだと遅いです。フィボナッチ数列の10万項までこの判定をして36秒でした。
どこが遅いか調べると、10^9の剰余の計算をするところです。これを避けるために、フィボナッチ数列の計算と並行して下9桁だけの計算をしましょう。これは32ビット整数の範囲で計算できるので高速です。

def gen_Fib_mod(d):
    a = 1
    b = 1
    while True:
        yield a
        a, b = b, (a + b) % d

これで7.7秒になりました。また、pandigitalなら9の倍数なので、これを判断に加えると、6.2秒になりました。
ここで上9桁の判定を加えると、

def is_upper_pandigital(n):
    a = digits(n)
    return is_pandigital(a[:9])

27秒になりました。たまにしか判定しないのに非常に遅いです。
桁数をフィボナッチ数列の公式を元にlogで計算して桁数を出して、10のべきで割ることにします。

def is_upper_pandigital(n, k):
    e = int((k * log((1 + sqrt(5)) / 2) - log(5)) / log(10) - 0.7)
    m = n / 10 ** (e - 9)
    if m >= 10 ** 9:
        m /= 10
    return is_pandigital(digits(m))

これで6.4秒でした。
これで全体の計算をすると65秒。


実数計算で概算して上9桁だけ取り出すと1.5秒でしたが、これだと誤差で上9桁の最後の数字が微妙になるとき、どうすれば正しくなるか分かりません。

def is_upper_pandigital(k):
    e = (k * log((1 + sqrt(5)) / 2) - log(5) / 2) / log(10)
    m = int(10 ** (e - int(e) + N))
    while m >= D:
        m /= 10
    return is_pandigital(digits(m))
from itertools import *
from math import sqrt, log

def digits(n):
    def gen_digits(n):
        while n:
            n, r = divmod(n, 10)
            yield r
    
    return list(gen_digits(n))

def gen_Fib():
    a = 1
    b = 1
    while True:
        yield a
        a, b = b, a + b

def gen_Fib_mod(d):
    a = 1
    b = 1
    while True:
        yield a
        a, b = b, (a + b) % d

def is_pandigital(ds):
    return sum(1 << d for d in ds) == 1022

def is_lower_pandigital(n):
    return n % 9 == M and is_pandigital(digits(n))

def is_upper_pandigital(n, k):
    e = int((k * log((1 + sqrt(5)) / 2) - log(5) / 2) / log(10) - 0.4)
    m = n / 10 ** (e - N)
    while m >= D:
        m /= 10
    return is_pandigital(digits(m))

def is_valid(x):
    f, f1 = x[1]
    return is_lower_pandigital(f1) and is_upper_pandigital(f, x[0])

N = 9
M = N * (N + 1) / 2 % 9
D = 10 ** N
print (x[0] for x in ifilter(lambda x: is_valid(x),
            izip(count(1), izip(gen_Fib(), gen_Fib_mod(D))))).next()