問題の概要
- d(n)をnの「真の約数」の和とする。
- 「真の約数」はnの約数のうちnを除外したものと定義する。
- d(a) == b かつ d(a) == bであるようなaとbをAmicable numbers(友愛数)と呼ぶ。
- 完全数(すなわちa == b)は除外する。
- 10000以下の友愛数の総和を求めよ。
解答
最初完全数を除外する条件を見逃していてハマった。
ちなみに問題文からはaが10000以下でbが10001以上の場合の処理方法が分からないが、とりあえず今回の範囲ではそのような友愛数のペアは出現しない。
何も考えずに解く
# constant ---- limit = 10000 # function ---- d <- function(n){ if(n == 1) return(0) if(n == 2) return(1) sum((1:(n-1))[n %% (1:(n-1)) == 0]) } # solve ---- f1 <- function(){ result <- 0 for(i in 2:(limit)){ d1 <- d(i) d2 <- d(d1) if(i == d2 && i != d1){ result <- result + i } } return(result) }
答えは出たがギリ1000msオーバーというところでお話にならない。
> microbenchmark::microbenchmark(f1()) Unit: milliseconds expr min lq mean median uq max neval f1() 779.8276 1026.256 1113.886 1091.904 1173.034 1619.54 100
約数の和の計算を高速化
d()
を少し修正してみる。
1つ約数が見つかった時点でもう1つ見つかっているという点に注意すれば約数の探索範囲はまでで良い。ただし平方数の場合の処理を忘れないようにする(書いてから気づいたけどforの外で処理すべきだ)。
d <- function(n){ result = 1 for(i in 2:sqrt(n)){ if(n %% i == 0){ result = result + i + n/i } if(i * i == n) { result = result - i } } return(result) }
これだけでもかなり早くなる。
> microbenchmark::microbenchmark(f1()) Unit: milliseconds expr min lq mean median uq max neval f1() 202.7227 245.5611 292.1344 279.2209 328.2863 478.2382 100
これはR3.4.0からデフォルトで有効になっているJITコンパイラの恩恵も大きく、試しにcompiler::enableJIT(0)
としてJITコンパイラをOFFにしてみると計算時間は倍になる。
> microbenchmark::microbenchmark(f1()) Unit: milliseconds expr min lq mean median uq max neval f1() 439.3626 525.4667 624.8876 600.091 701.4945 1021.952 100
Rcppの力を借りる
多少早くなったといってもRのforなのでまだまだ遅い。試しに一部をRcppで置き換えてみる。
#include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] int dC(int n) { int result = 1; if (n <= 1) { return 0; } else if (n == 2) { return 1; } else { for(int i = 2; i*i <= n; i++){ if(n % i == 0){ result += i + n/i; } if(i*i == n){ result -= i; } } } return result; }
> microbenchmark::microbenchmark(f2()) Unit: milliseconds expr min lq mean median uq max neval f2() 34.38936 39.30994 44.53533 42.54718 44.37687 117.9409 100
大分早い。いっそもう全部置き換えると…
int LIMIT = 10000; // [[Rcpp::export]] int getans() { int result = 0; for(int i = 2; i <= LIMIT; i++){ int d1 = dC(i); int d2 = dC(d1); if(i == d2 & i != d1){ result += i; } } return result; }
> microbenchmark::microbenchmark(getans()) Unit: milliseconds expr min lq mean median uq max neval getans() 4.410787 4.716603 5.319099 5.289165 5.732805 7.585598 100
最初からC++でやれという感じに。