Problem 12 三角数の約数　解答編

偉そうに解答編と銘打ってみます。
id:selvaggio:20080328の続きです。

素直に

まずCでの素直なコード。

#include <stdio.h>

int divs(int n){
  int i=1, ret = 0;
  for(; i<=n; ++i)
    if(n%i==0) ret++;
  return ret;
}

int tri(int n){
  return n*(n+1)/2;
}

int main()
{
  int i=1,t;
  for(;;i++){
    t=divs(tri(i));
    if(t>500){
      printf("%d %d %d\n",i, tri(i), t);
      break;
    }
  }
}

divsが約数をカウントする関数。
triが三角数を計算する関数。
mainは、三角数の約数をカウントし、500を超えたらそれを表示してbreakする。

このプログラムを実際に走らせると分かるけど、答えはなかなかでない。
divsが遅いからだ。

divsは何故遅いのか？

例えば42の約数を考えてみる。約数は1,2,3,6,7,14,21,42の8つだ。
この約数の列をじっと見つめると、ある性質に気づく。
後半に行くに従って、約数が疎らになっている。
視覚的に分かりやすく表現するとこんな感じ。
123..67......14.......21..................42
計算時間の大部分はこの後半の"."に割かれていることになる。

改善策

約数の列をもう一度じっと見つめると、もう一つの性質に気づく。
こんな具合に。
42(元の数) = 1 * 42 = 2 * 21 = 3 * 14 = 6 * 7 = 7 * 6 = 14 * 3 = 21 * 2 = 42 * 1(二つの約数の積)

これは偶然ではない。
「 $d$ が $N$ の約数のとき、 $e=N/d$ も約数である」という性質があるからだ。
この性質を利用すると、divsは劇的に効率が良くなる。

どうするかといえば、
上の式で言えば、6*7のところまでで全ての約数が式に登場していて、後半は鏡写しになっているのを利用する。
鏡写しになる条件は $N=d*e$ の $d$ と $e$ が $d>e$ のときだ。

この不等式を解くと、 $d>sqrt{N}$ となる。
つまり、ループは $d \le sqrt{N}$ の間まわせばいい。
ただし、一つ約数を見つけたときはペアとなる約数があるので、カウントは2つずつ増える。
また、36の約数の6のように、自分自身がペアになる場合があるので注意。

コード

新しいdivsはこうなる。

int divs(int n){
  int i=1, ret = 0;
  for(; i<sqrt(n); ++i)
    if(n%i==0) ret+=2;
  if(i*i==n) ret++;
  return ret;
}

for文でペアを数え、その後のifで自分自身とペアになっている約数をカウントしている。
これだけでもdivsは相当速くなる。
オーダーで言えば、 $O(N) \to O(sqrt{N})$ だ。劇的。

まだ遅い？

ここで満足しそうだけど、まだ改善の余地はある。
でも長くなってきたのでまた次回ということで。