Problem 95

R Project Euler

いずれの要素も1,000,000を超えない社交数(c.f. 社交数 - Wikipedia．ちなみに答が書いてある．)のうち，要素数が最大であるものの最小の要素を答えろ，という問題．
計算時間を短縮する方法がさっぱり思い付かない．480秒もかかった．またそのうちやろう(≒多分もうやらない)．
NGワード: ていうかRじゃなくね？，もう全部Cでいいじゃん
inlineの使い方がほんのり分かったのは収穫．ほんのりしか分かってないので難しいことはできないし説明もできない．というかそもそもCが読めない書けない．
ちなみにRだけでやると1時間半近くかかる．約数の和を計算する部分を置き換えただけで10倍くらい早くなっているということ．やっぱProject EulerはRでやるもんじゃないなと思いました．

library(inline)
code <- "
  int n = INTEGER(x)[0];
  SEXP sum;
  PROTECT(sum = NEW_NUMERIC(1));
  REAL(sum)[0] = 1;
  if(n == 1) return(sum);
  int r = (int)sqrt(n);
  int f, step;
  if(r * r == n){
    REAL(sum)[0] = 1 + r;
    r -= 1;
  } else {
    REAL(sum)[0] = 1;
  }
  if(n % 2 == 1){
    f = 3;
    step = 2;
  } else {
    f = 2;
    step = 1;
  }
  while(f <= r){
    if( n % f == 0 ){
      REAL(sum)[0] += f + n/f;
    }
    f += step;
  }
  UNPROTECT(1);
  return(sum);
"
sum.of.div <- cfunction(signature(x = "numeric"), code)

limit <- 1e6
check.list <- logical(limit)
max.chain <- numeric(0)
for(i in 3:limit){
  if(check.list[i]) next
  n <- i
  temp.chain <- n
  repeat{
    check.list[n] <- TRUE
    n <- sum.of.div(as.integer(n))
    temp.chain <- c(temp.chain, n)
    if(n==1 || n > limit) break
    if(check.list[n]){
      temp.chain <- temp.chain[((1:length(temp.chain))[temp.chain==n][1]):length(temp.chain)]
      if(length(temp.chain) > length(max.chain)){
        max.chain <- temp.chain
      }
      break
    }
  }
}
min(max.chain)