Projecet Euler 12 割り算をせずに約数の個数を求める。

問題

三角数の数列は自然数の和で表わされ, 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に変えれば約数の個数が求められるので、それを利用するともっと早くなると思う。