Project Euler Problem 012 - つまみ食う

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

The sequence of triangle numbers is generated by adding the natural numbers.
So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

  1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Let us list the factors of the first seven triangle numbers:

   1: 1
   3: 1,3
   6: 1,2,3,6
  10: 1,2,5,10
  15: 1,3,5,15
  21: 1,3,7,21
  28: 1,2,4,7,14,28
  
We can see that 28 is the first triangle number to have over five divisors.
What is the value of the first triangle number to have over five hundred divisors?

約数が500個を超える最初の三角数を求める。

まず約数の数を数えるコードが必要。単純に割れる数をカウントする。こんな感じ。

import math

def divisorlen(n):
    if n == 1: return 1
    
    m = math.sqrt(n)
    result = 2
    
    i = 2
    while i <= m:
        if n % i == 0:
            result += 2
        i += 1
    else:
        if i - 1 == m:
            result -= 1
    
    return result

あとは閾値を超えるまで三角数の約数の数を求めていく。

import datetime
import itertools    

def triangle(n):
    return n * (n + 1) / 2

def euler012():
    begin = datetime.datetime.now()
    
    for i in itertools.count(1):
        x = triangle(i)
        l = divisorlen(x)
        if l > 500:
            print 'f(%d) = %d, len=%d' % (i, x, l)
            break
        
    end = datetime.datetime.now()
    print end - begin

しかしこのコードは相当遅い。14秒くらいかかる。

改善

さっきの方法だと大きい数の約数の数を求める必要があって、それがボトルネックとなる。少し工夫すれば比較的小さい数を使って同じ計算ができる。とりあえず説明を書くけど、数学ができないのできちんと説明できないと、さきに言い訳をしておく。

三角数の式 n(n+1)/2 に注目する。

n と n+1 は約数が重複しない（これを「互いに素」と言うらしい）
互いに素な二つの数字は約数が重複しないので

  n(n+1) の約数の数 = nの約数の数 * n+1の約数の数

さらに

n と n+1 のどちらかは偶数である。
偶数aの約数の集合をAとすると、a/2の約数の集合はAのサブセットとなる（なので互いに素な状態は保たれる）

なので、n(n+1)/2 の約数の数を求めるには

  nが偶数の場合: n/2の約数の数 * n+1の約数の数
  nが奇数の場合: nの約数の数 * (n+1)/2の約数の数

となる。

この記事が同じような内容を分かりやすく説明していた。
http://d.hatena.ne.jp/keisukefukuda/20100729/p1
約数の数を求めるのにキャッシュを使ったり、僕がこのあとに書くコードよりもさらに工夫されている。なるほど！！！！

そのことを踏まえてコード化すると

import datetime
import itertools    

def triangle(n):
    return n * (n + 1) / 2

def euler012():
    begin = datetime.datetime.now()
    
    a, b = 0, 1
    for i in itertools.count(1):
        
        if i % 2 == 0:
            a = divisorlen(i/2)
            b = divisorlen(i+1)
        else:
            a = divisorlen(i)
            b = divisorlen((i+1)/2)
            
        if a * b > 500:
            x = triangle(i)
            print 'f(%d) = %d, len=%d' % (i, x, a * b)
            break
        
    end = datetime.datetime.now()
    print end - begin

だいぶん早くなる。0.4秒くらい。

でも、まだ改善の余地がある。n+1 は次の三角数の n なので計算量は半分になる。

import datetime
import itertools    

def triangle(n):
    return n * (n + 1) / 2

def euler012():
    begin = datetime.datetime.now()
    
    a, b = 0, 1
    for i in itertools.count(1):
        
        a = b
        if i % 2 == 0:
            b = divisorlen(i+1)
        else:
            b = divisorlen((i+1)/2)
            
        if a * b > 500:
            x = triangle(i)
            print 'f(%d) = %d, len=%d' % (i, x, a * b)
            break
        
    end = datetime.datetime.now()
    print end - begin

答え: f(12375) = 76576500, len=576
実行時間: 0.219535秒くらい

f(1)から始めるのではなく途中から始めた方がいいかもしれないけど、じゃあ何番から始めるのが妥当なのかの判断が出来なかった。