Problem 150

Problem 150 - Project Euler

正負の整数からなる三角配列から、含まれる数の和が最小となる部分三角配列を探したい。
例を示そう(上記リンク先原文中図参照)。例では枠で括った部分の配列の総和が-42であり条件を満たすことが容易に確認できる。
1000行の三角配列を作りたい。そこで、±2^19の範囲の擬似乱数を500500個、次の乱数生成法(線形合同法)により生成する。

t := 0
for k = 1 up to k = 500500:
  t := (615949*t + 797807) modulo 2**20
  s[k] := t - 2**19

これにより生成される乱数は、s[1] = 273519、s[2] = -153582、s[3] = 450905、...etcとなる。
三角配列には次のように数字を入れていく。

        s[1]
      s[2] s[3]
   s[4] s[5] s[6]
s[7] s[8] s[9] s[10]

配列のある要素を始点とする部分三角配列は、下に向かって広がっていくものと考える(つまり、2番目の行は2つの要素を含み、その次の行は3つの要素を含み、というように続けて部分三角配列を作る)。
「部分三角配列の和」は部分三角配列に含まれる要素の総和であるとする。
最小の部分三角配列の和を求めよ。

Problem 149と似たような感じでやった。

t <- 0
s <- numeric(500500)
for(k in 1:500500){
  t <- (615949 * t + 797807) %% (2^20)
  s[k] <- t - (2^19)
}
s[1:3]

tri <- matrix(TRUE, nrow=1000, ncol=1000)
tri[lower.tri(tri)] <- FALSE
tri[tri] <- s
tri <- t(tri)

## 累積和に変更
tri.cs <- tri
for(i in 1:1000){
  for(j in 1:i){
    tri.cs[i, j] <- sum(tri[i, 1:j])
  }
}

smin <- 0
for(i in 1:1000){
  for(j in 1:i){  ## 始点: [i, j]
    s <- tri.cs[i, j] ## s: 累積和
    if(j > 1) s <- s - tri.cs[i, j-1]
    l <- 1 ## 始点からの行数
    for(k in (i+1):1000){
      s <- s + tri.cs[k, j+l]
      if(j > 1) s <- s - tri.cs[k, j-1]
      if(s < smin) smin <- s
      l <- l + 1
    }
  }
}
smin

最後まで回してないけどたぶん1時間以上かかる。
が、Rcppの力を借りると20秒足らず。ううむ…。
CやC++で書くより楽(というか僕はCやC++書けない)だしRで解いたという扱いでいいか…。

library(Rcpp)
library(inline)

make.s <- cxxfunction(body="
IntegerVector s(500500);
int t = 0;
for(int k = 0; k < 500500; k++){
  t = (615949 * t + 797807) % (int)pow(2, 20);
  s[k] = t - (int)pow(2, 19);
}
return wrap(s);
", plugin = "Rcpp")

s <- make.s()
s[1:3]

tri <- matrix(TRUE, nrow = 1000, ncol = 1000)
tri[lower.tri(tri)] <- FALSE
tri[tri] <- s
tri <- t(tri)

## 累積和に変更
tri.cs <- tri
for(i in 1:1000){
  for(j in 1:i){
    tri.cs[i, j] <- sum(tri[i, 1:j])
  }
}

find.ans <- cxxfunction(signature(t = "integer"),"
IntegerMatrix tri(t);
int smin = 0;
for(int i=1; i<=1000; i++){
  for(int j=1; j<=i; j++){
    int s = tri(i, j);
    if(j > 1) s -= tri(i, j-1);
    int l = 1;
    for(int k=(i+1); k<=1000; k++){
      s += tri(k, j+l);
      if(j > 1) s -= tri(k, j-1);
      if(s < smin) smin = s;
      l += 1;
    }
  }
}
return wrap(smin);
", plugin = "Rcpp")

find.ans(tri.cs)