オイラー関数のメモ化 素数表つき

オイラー関数のメモ化

まずはメモ化しない愚直なコード

# euler.py

def gcd(x, y):
    while y != 0:
        x,y = y,x%y
    return x

def euler(n):
    if n==1: return 1
    cnt = 0
    for i in xrange(1, n):
        if gcd(n, i) == 1:
            cnt += 1
    return cnt

def main():
    for n in xrange(1, 10000):
        euler(n)

if __name__ == '__main__':
    main()

実行：

$ time python euler.py 

real	1m45.373s
user	1m41.762s
sys	0m0.080s

Psycoの"import psyco; psyco.profile()"を使って高速化すると：

$ time python euler.py 

real	0m7.349s
user	0m7.000s
sys	0m0.004s

次にちょっとかしこいメモ化なしのコード

# euler2.py
def gcd(x, y):
    while y != 0:
        x,y = y,x%y
    return x

def euler(n):
    if n==1: return 1
    cnt = 0
    for i in xrange(1, n):
        g = gcd(n, i)
        if g == 1:
            cnt += 1
        else:
            if gcd(n/g, i) == 1:
                return euler(n/g) * euler(i)
    return cnt

def main():
    for n in xrange(1, 10000):
        euler(n)

if __name__ == '__main__':
    main()

実行：

$ time python euler2.py 

real	0m21.795s
user	0m21.221s
sys	0m0.028s

Psyco付きだと：

$ time python euler2.py 

real	0m1.543s
user	0m1.516s
sys	0m0.000s

メモ化コード

# euler_memoize.py
# euler(p) = p-1 (p is prime)
# euler(m*n) = euler(m) * euler(n) (m,nは互いに素)

def gcd(x, y):
    while y != 0:
        x,y = y, x%y
    return x

# 樹状再帰するメモ化オイラー関数
cache = {0:0, 1:1, 2:1}
def euler(n):
    global cache
    if n in cache:
        return cache[n]
    cnt = 0
    for i in xrange(1, n):
        g = gcd(i, n)
        if g == 1:
            cnt += 1
        else:
            if gcd(n/g, i) == 1:
                #print "phi(%d) = phi(%d)*phi(%d)" %(n,n/g,i)
                cnt = euler(n/g) * euler(i)
                cache[n] = cnt
                return cache[n]
    cache[n] = cnt
    return cache[n]
    
def main():    
    for x in xrange(1, 10000):
        euler(x)

if __name__ == '__main__':
    main()

実行：

$ time python euler_memoize.py 

real	0m14.088s
user	0m13.897s
sys	0m0.008s

Psyco付き：

$ time python euler_memoize.py 

real	0m0.998s
user	0m0.960s
sys	0m0.008s

速くなった。最後に、明らかに、引数が素数の場合に時間がかかる。オイラー関数(p) = p-1 (p:素数) なので、予めメモしておけばさらに早くなるはず。

# euler(p) = p-1 (p is prime)
# euler(m*n) = euler(m) * euler(n) (m,nは互いに素)

#import psyco; psyco.profile()

def gcd(x, y):
    while y != 0:
        x,y = y, x%y
    return x

# 樹状再帰するメモ化オイラー関数
cache = {0:0, 1:1, 2:1}
def euler(n):
    global cache
    if n in cache:
        return cache[n]
    cnt = 0
    for i in xrange(1, n):
        g = gcd(i, n)
        if g == 1:
            cnt += 1
        else:
            if gcd(n/g, i) == 1:
                #print "phi(%d) = phi(%d)*phi(%d)" %(n,n/g,i)
                cnt = euler(n/g) * euler(i)
                cache[n] = cnt
                return cache[n]
    cache[n] = cnt
    return cache[n]

def fill_prime():
    global cache
    table = range(10000)
    table[1] = 0
    for i in xrange(2, 100):
        for k in xrange(i*2, 10000, i):
            table[k] = 0
    for p in [x for x in table if x != 0]:
        cache[p] = p-1
    
def main():
    fill_prime()
    
    for x in xrange(1, 10000):
        euler(x)

if __name__ == '__main__':
    main()

実行：

$ time python euler_memoize.py 

real	0m0.445s
user	0m0.440s
sys	0m0.008s

Psycoつき

$ time python euler_memoize.py 

real	0m0.113s
user	0m0.112s
sys	0m0.004s

更に速くなった。