正負の整数からなる三角配列から、含まれる数の和が最小となる部分三角配列を探したい。
例を示そう(上記リンク先原文中図参照)。例では枠で括った部分の配列の総和が-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)