下の表において、縦、横、あるいは斜め方向に連続して(注:長さは問わず、途中で曲がってはいけない)隣接する数の和で、最大となるものは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使おう。