Project Euler Problem 21

R Project Euler

問題の概要

Problem 21 - PukiWiki

  • 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つ見つかっているという点に注意すれば約数の探索範囲は\sqrt{n}までで良い。ただし平方数の場合の処理を忘れないようにする(書いてから気づいたけど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++でやれという感じに。

参考