問題
三角数の数列は自然数の和で表わされ, 7番目の三角数は 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28 である. 三角数の最初の10項は:
1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...
となる.
最初の7項について, その約数を列挙すると, 以下のとおり.
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
これから, 7番目の三角数である28は, 5個より多く約数をもつ最初の三角数であることが分かる.
では, 500個より多く約数をもつ最初の三角数はいくつか.
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2012
数学的考察
三角数というと難しげだが、要は三角数とは1からnまでの自然数の和で表される数である。
1からnまでの自然数の和は、以下の公式で表される。
n * (n+1) /2
そして、nと(n+1)とは互いに素（共通の素因数を持たない。）であるから、n/2とn+1またはnと(n+1)/2も互いに素である。n/2とn+1とすべきか、nと(n+1)/2とすべきかはいずれが偶数かによる。
互いに素な自然数n1,n2をかけあわせた数の約数の個数は、n1の約数の個数とn2の約数の個数とを掛け合わせた数に等しい。
回答方針1
上記のような数学的考察に基づき、以下のようなアルゴリズムで本問の回答を行う。
1. mymathを使って素数のリストを作る。
mymath.py
http://qiita.com/cof/items/45d3823c3d71e7e22920
2. 約数の個数が500以内であるうちはループ ループ内では変数nを1 から順に増加させる。
3. n+1について偶数か否かを判定した上で、素数リストを使ってn+1または(n+1)/2の素因数分解を行う。
4. 得られた素因数を下にn+1の約数の数を求め、以前に記憶していたnの約数の数と掛け合わせる。
5. かけあわせた数を約数の数として記憶、2に戻る。
コード
import mymath
def cof():
pri = mymath.get_primes(10**4)
n = 1
num_fac = 1
lim = 500
n1_fac = mymath.get_number_of_factor(n,pri)
while num_fac < lim:
n += 1
if (n+1)%2:
n2_fac = mymath.get_number_of_factor(n+1,pri)
else:
n2_fac = mymath.get_number_of_factor((n+1)/2,pri)
num_fac = n1_fac * n2_fac
n1_fac = n2_fac
print n * (n+1) / 2
問題の発覚と新たな気づき
実行してみたところ、1回の実行に1.62秒と、うんこレベルで遅い。
どうにかして高速化できないかと考えてみたところ、そもそも因数分解をしていることが遅い原因ではないか、ということに気がついた。
エラトステネスと同様の手法によって、割り算をせずに、各数の約数を求められることに気がついた。
回答方針2
より詳しくは、
1. それぞれの数ごとにリストを作る。
[[],[],・・・]
2. 1から適当な数までループし、対象の数の倍数の位置に当該対象の数を要素として追加する。
例えば、2が対象の数であれば、
[[].[1],[1],[1],[1]←倍数の4に対応するリストに2を追加する、以下、6,8,10と同様。
3. 2を繰り返すと、あら不思議、対応する数の真の約数リストが出来上がります。
[[].[1],[1],[1],[1,2],[1],[1,2,3],[1],[1,2,4],[1,3],[1,2,5],...]
4. この約数リストを基に、回答方針1にあるnの約数の個数等を求める。
コード
以下の関数をmymathに追加。
def factor_seq(max):
ret = [[0]]
for i in range(max):
ret += [[1]]
seq = range(max+1)
for i in seq[2:max//2+1]:
for j in seq[i*2::i]:
ret[j] += [i]
return ret
上記を使ってコードディング
def cof2():
fq = mymath2.factor_seq(10**5)
n = 1
num_fac = 1
lim = 500
n1_fac = len(fq[n])
while num_fac < lim:
n += 1
if (n+1)%2:
n2_fac = len(fq[n+1])
else:
n2_fac = len(fq[(n+1)/2])
num_fac = n1_fac * n2_fac
n1_fac = n2_fac
print n * (n+1) / 2
結果、640msと、一応なんとか我慢できるレベルにはなった。たぶん、factor_seqでret[j] += [i]をret[j] += 1に変えれば約数の個数が求められるので、それを利用するともっと早くなると思う。