2 〜 n までの整数の約数の個数

2 ~ n までの整数を素因数分解したとき，約数（素数） p が何個あるか数える。

p = 7 のとき 35 = 5 * 7 なので 1 個，49 = 7 * 7 なので 2 個である。

素朴にプログラムしてみると，以下のようになる（素朴すぎるということがすぐわかる）。

> count_factors0 = function(n, p) {
+   count = 0
+   for (i in 2:n) {
+     while (i %% p == 0) {
+       count = count + 1
+       i = i %/% p
+     }
+   }
+   return(count)
+ }

2 〜 50 の約数 7 の個数は，

> count_factors0(50, 7)
[1] 8

ただし，大きな n に対して求めようとすると，とんでもない時間が掛かる。

> system.time({
+   print(count_factors0(50000000, 7))
+ })
[1] 8333328
   ユーザ   システム       経過  
    14.356      0.030     14.386 

for ループと while ループの使い方が悪いのだ。以下のようにすれば瞬殺である。

> count_factors2 = function(n, p) {
+   count = 0
+   while (n > 0) {
+     count = count + n %/% p
+     n = n %/% p
+   }
+   return(count)
+ }

> system.time({
+   print(count_factors2(50000000, 7))
+ })
[1] 8333328
   ユーザ   システム       経過  
     0.004      0.000      0.004 

ちょっとわかりにくいが，以下のように再帰関数で定義するとわかりやすいかな。

> count_factors1 = function(n, p) {
+   if (n == 0) {
+     return(0)
+   } else {
+     return(n %/% p + Recall(n %/% p, p))
+   }
+ }

> system.time({
+   print(count_factors1(50000000, 7))
+ })
[1] 8333328
   ユーザ   システム       経過  
     0.005      0.000      0.005 

外積の和

外積は R では outer(a, b, 演算子) で計算できる。

c(1, 2, 3, 4, 5) と c(2, 3, 4) で，演算子を "*" （積）とすると，以下のようになる。

> x = outer(1:5, 2:4, "*")
> x
     [,1] [,2] [,3]
[1,]    2    3    4
[2,]    4    6    8
[3,]    6    9   12
[4,]    8   12   16
[5,]   10   15   20

結果の行列の各要素の和は sum()  で求められる。

> print(sum(x))
[1] 135

さて，ここで，この結果を得るための最良の方法を探索しよう。

長さが共に n = 1e4 の正規乱数ベクトル a, b の sum(outer(a, b, "*")) を求めることとする。

outer() を使えば以下のようになる。答えは 2159.779 であり，0.309 秒ほどかかった。これでも十分速い。

> system.time({
+   x = outer(a, b, "*")
+   print(sum(x)) # 2159.779
+ })
[1] 2159.779
   ユーザ   システム       経過  
     0.309      0.150      0.459 

outer() がやっていることを R で素直（？）に書けば以下のようになる。

> system.time({
+   sum.of.outer.multiple = 0
+   for (i in seq_along(a)) {
+     for (j in seq_along(b)) {
+       sum.of.outer.multiple = sum.of.outer.multiple + a[i] * b[j]
+     }
+   }
+   print(sum.of.outer.multiple)
+ })
[1] 2159.779
   ユーザ   システム       経過  
     6.773      0.019      6.818

20  倍も遅い。

内側の for ループは，a[i] * (b[1] + b[2] + ... + b[n]) を計算している訳だが，ループでは (b[1] + b[2] + ... + b[n]) はいつも同じなので，前もって計算しておいてそれを使えば，ループが省けるし計算量も削減できる。

> system.time({
+   sum.b = sum(b)
+   sum.of.outer.multiple = 0
+   for (a.i in a) {
+     sum.of.outer.multiple = sum.of.outer.multiple + a.i * sum.b
+   }
+   print(sum.of.outer.multiple)
+ })
[1] 2159.779
   ユーザ   システム       経過  
     0.003      0.000      0.003 

いやいや，待って待って。この for ループも (a[1] + a[2] + ... + a[n]) * sum.b を計算していますね。

な〜んだ。sum.a = a[1] + a[2] + ... + a[n] としておいて，sum.a * sum.b で答えが出ますね。

> system.time({
+   print(sum(a)*sum(b))
+ }) # 0.000
[1] 2159.779
   ユーザ   システム       経過  
         0          0          0 

瞬殺でした。outer() を使う必要はなかった。