Problem 149

Problem 149 - Project Euler

下の表において、縦、横、あるいは斜め方向に連続して(注:長さは問わず、途中で曲がってはいけない)隣接する数の和で、最大となるものは16(=8 + 7 + 1)であるということは容易に確認できる。

-2 5 3 2
9 -6 5 1
3 2 7 3
-1 8 -4 8

より大きなスケールの表に関して、このような和を探索しよう。
まず、400万個の擬似乱数を、ラグ付きフィボナッチ法と呼ばれる方法で生成しよう。
For 1 ≦ k ≦ 55, .
For 56 ≦ k ≦ 4000000, .
ここで、となる。
次に、2000×2000の表を容易する。そして今作った乱数の最初の2000個を最初の行に、次の2000個を次の行に、といった具合に数字を入れていって表を作成する。
最後に、作成した表の中から縦、横、斜めに連続して隣接する数の和で最大のものを見つけよ。

2000行と2000列、それと斜め方向に3999列×2の各々の数列について最大部分列和を求め、その中で最大のものを見つければいい。
斜め方向が若干面倒だが、左上半分だけ探索する関数を作って、表のほうを回転させることで対応した。

s <- numeric(4e6)
k <- 1:55
s[k] <- ((100003 - 200003*k + 300007*k^3) %% 1000000) - 500000
for(k in 56:4e6){
  s[k] <- ((s[k-24] + s[k-55] + 1000000) %% 1000000) - 500000
}

smat <- matrix(s, nrow=2000)

mss <- function(a){
  s <- mss <- 0
  for(i in 1:length(a)){
    s <- s + a[i]
    if(s < 0) s <- 0
    if(mss < s) mss <- s
  }
  mss
}

limit <- 2000
m <- 0
for(i in 1:limit){ ## 縦横
  m <- max(m, mss(smat[i,]), mss(smat[,i]))
}

## 斜め(/)の上半分を調べる
slash <- function(mat){
  m <- 0
  for(i in 2:nrow(mat)){
    s <- mt <- 0
    for(j in 1:(i-1)){
      s <- s + mat[i-j, j]
      if(s < 0) s <- 0
      if(mt < s) mt <- s
    }
    if(m < mt) m <- mt
  }
  m
}
m <- max(m, slash(smat),
         slash(smat[,2000:1]),
         slash(smat[2000:1,]),
         slash(smat[2000:1, 2000:1])
         )

上記コードで大体2分かかる。そもそも乱数表を生成している時点で数十秒かかっており、なんかどうしようもない感じがする。
もっと上手いやり方やRにあったコーディングの仕方を工夫すればいいのかもしれないが、とりあえずRcppパッケージとinlineパッケージを使ってみることにした。

library(Rcpp)
library(inline)

Lagged.Fibonacci.Generator <-
  cxxfunction(body ="
IntegerVector s(4000000);
for(int k=1; k<=55; k++){
  s[k-1] = ((100003L - 200003L*k + 300007L*k*k*k) % 1000000L) - 500000L;
}
for(int k=55; k<4000000; k++){
  s[k] = ((s[k-24] + s[k-55] + 1000000L) % 1000000L) - 500000L;
}
return wrap(s);
", plugin = "Rcpp")
s <- Lagged.Fibonacci.Generator()
smat <- matrix(s, nrow=2000)

## 縦横
mss <- cxxfunction(signature(a = "integer"),"
IntegerVector aa(a);
int s = 0;
int mss = 0;
int n = aa.size();
for(int i=0; i<n; i++){
  s += aa[i];
  if(s < 0) s = 0;
  if(mss < s) mss = s;
}
return wrap(mss);
", plugin = "Rcpp")

## 斜め(/)の上半分
slash <- cxxfunction(signature(a = "integer"),"
IntegerMatrix mat(a);
int m = 0;
for(int i=0; i<2000; i++){
  int s = 0; int mt = 0;
  for(int j=0; j<=i; j++){
    s += mat(i-j, j);
    if(s < 0) s = 0;
    if(mt < s) mt = s;
  }
  if(m < mt) m = mt;
 }
return wrap(m);
", plugin = "Rcpp")

m <- 0
for(i in 1:limit){ ## 縦横
  m <- max(m, mss(smat[i,]), mss(smat[,i]))
}
max(m, slash(smat),
    slash(smat[,2000:1]),
    slash(smat[2000:1,]),
    slash(smat[2000:1, 2000:1])
    )

全部実行して、つまりC++部分のコンパイル時間込みで11秒程度、コンパイル済みの関数を使って計算するだけだと500msちょい上回る程度で答えが出る。
まあ問題はほとんどRのコードじゃないってところなんだけど…。

結論

Rcpp使おう。