> 「データ数」ってのはちょっと変で（データは複数形だから）「サンプル数」ってのが正しいのかと思ってた．なんでいけないのかね．自分は「データの数がさあ」とか普通に怒鳴ってるけど．統計学者というのはよくわからないことにこだわる人たちではある．
科学者（専門家）が正しい用語を使うこと，および他者に対して正しい用語を使うように要請することを，「統計学者というのはよくわからないことにこだわる人たちではある」と茶化すのはどうかなあ。
pi を，試行回数 n = 10, 100, 1000, ..., 1000000 として，モンテカルロ法で求めよという単純な問題。
ずっと前にここにも書いたが，R では基本的に
mean(1 >= colSums(matrix(runif(2*n)^2, 2)))*4
の 1 行で書ける（複数の関数を使うが）
junk = sapply(10^(1:6), function(n) {
pi0 = mean(1 >= colSums(matrix(runif(2*n)^2, 2)))*4
cat(sprintf("n = %7i pi = %.9f error = %10.7f\n", n, pi0, pi-pi0))
}
n = 10 pi = 3.200000000 error = -0.0584073
n = 100 pi = 3.120000000 error = 0.0215927
n = 1000 pi = 3.252000000 error = -0.1104073
n = 10000 pi = 3.150000000 error = -0.0084073
n = 100000 pi = 3.142480000 error = -0.0008873
n = 1000000 pi = 3.144192000 error = -0.0025993
なお，
「ランダムな点を増やしていくと、段々と Πに近づいていくはずです」
というのは，一見正しそうであるが，一時的に遠ざかることもあるということで，実行例のような n が小さいときの値の方が誤差が少ないということもある（あたりまえなんだ）。
また，
「精度は単精度で構いません」
とは，時代遅れ。今時，単精度でプログラムすることなんかないだろう。
ガラケーを使っている人が出した暗号というのが，ヒントなんだろうとピンときた。
2 タッチ入力かと思ったが，ちょと違った。
ケータイ入力の方だ。
まずは，平文化
decoding = function(x) {
char = list(c("あ", "い", "う", "え", "お", "ぁ", "ぃ", "ぅ", "ぇ", "ぉ"),
c("か", "き", "く", "け", "こ"),
c("さ", "し", "す", "せ", "そ"),
c("た", "ち", "つ", "て", "と", "っ"),
c("な", "に", "ぬ", "ね", "の"),
c("は", "ひ", "ふ", "へ", "ほ"),
c("ま", "み", "む", "め", "も"),
c("や", "ゆ", "よ", "ゃ", "ゅ", "ょ"),
c("ら", "り", "る", "れ", "ろ"),
c("わ", "を", "ん", "ゎ"))
res = NULL
result = sapply(x, function(z) {
z1 = z%/%10
if (z1 == 0) z1 = 10
z2 = z%%10
res = c(res, char[[z1]][z2])
})
return(res)
}
decoding(c(11,32,41,25,46,35,92,33,71,65,52,43,12,44,15,32,14,44,23,94))
decoding(c(25,94,21,91,75,94,03,91,23,61,25,55,65,13,65,13))
暗号化も考えれば，文字と数字の対応ベクトルを作っておけばよい。
2 つの関数に同じデータを書くのは気が進まないので，暗号化か平文化かを引数で指定するようにしよう。
なお，作者の暗号法では，「ぉ」をどのように表すかが不定であることが明らかになったが，「ぉ」は 10 で表すようにしておく。
結果は関数の中で cat によって書き出す仕様にした。
ango = function(x, way=c("decoding", "encoding")) {
char = c("あ", "い", "う", "え", "お", "ぁ", "ぃ", "ぅ", "ぇ", "ぉ",
"か", "き", "く", "け", "こ",
"さ", "し", "す", "せ", "そ",
"た", "ち", "つ", "て", "と", "っ",
"な", "に", "ぬ", "ね", "の",
"は", "ひ", "ふ", "へ", "ほ",
"ま", "み", "む", "め", "も",
"や", "ゆ", "よ", "ゃ", "ゅ", "ょ",
"ら", "り", "る", "れ", "ろ",
"わ", "を", "ん", "ゎ")
num = c(11:19, 10, 21:25, 31:35, 41:46, 51:55, 61:65, 71:75, 81:86, 91:95, "01", "02", "03", "04")
way = match.arg(way)
if (way == "decoding") {
x = as.character(x)
result = sapply(x, function(z) {
if (nchar(z) == 1) z = sprintf("0%s", z)
char[grep(z, num)]
})
cat(paste(result, collapse=""))
} else {
x = unlist(strsplit(x, ""))
result = sapply(x, function(z) {
sprintf("%02i", as.integer(num[grep(z, char)]))
})
cat(result)
}
}
ango(c(11,32,41,25,46,35,92,33,71,65,52,43,12,44,15,32,14,44,23,94))
ango(c(25,94,21,91,75,94,03,91,23,61,25,55,65,13,65,13))
ango("りょうかいしました", way="encoding")
ango("あしたのおひるあいてます", way="e")
ango(c(92, 86, 13, 21, 12, 32, 71, 32, 41), way="d")
ango(c(11, 32, 41, 55, 15, 62, 93, 11, 12, 44, 71, 33))
待ち合わせの時間と場所を表しているという
ヒント：2x+3
解答
2x+3=169 で，x=83=0x53 ASCII コードだろう
R なら簡単
intToUtf8 関数だけで答が出る
x = c(169, 147, 161, 173, 145, 161, 67, 159, 149, 67, 171, 133, 153, 133, 161, 169, 133, 159, 67, 137, 147, 161, 173, 151, 161, 173, 67, 139, 141)
intToUtf8((x-3)/2)
0 ～ 9 までの数字をすべて使った割り算の筆算のうち，「割られる数」が最も小さくなるような場合を求め，その式を答えよ。
なお、「割られる数」と「割る数」はともに正の整数とし，整数の商と余りを求めるものとする。
calc = function(numerator, denominator) {
func = function(n) {
return(as.numeric(unlist(strsplit(as.character(n), ""))))
}
numerator.a = func(numerator)
denominator.a = func(denominator)
ans = numerator%/%denominator
ans.length = nchar(as.character(ans))
numerator.length = nchar(numerator)
i = 1
j = 0
num = numerator.a[i]
prod = num.list = vector("list", ans.length)
for (k in seq_len(ans.length)) {
if (num == 0)
break
while (denominator > num) {
i = i + 1
if (i > numerator.length)
break
num = 10 * num + numerator.a[i]
}
j = j + 1
num.list[[j]] = func(num)
ans[j] = num%/%denominator
prod[[k]] = func(denominator * ans[j])
num = num%%denominator
}
ans = c(ans, denominator.a, func(num), unlist(prod), unlist(num.list))
return(length(table(ans)) == 10)
}
found = FALSE
for (i in 1:1000) {
for (j in 1:i) {
if (calc(i, j)) {
cat(i, "/", j, "\n")
found = TRUE
break
}
}
if (found)
break
}
190 ÷ 78 かな。
各人(name)の拘束日時（date, begin, end) について以下のような予定表がある（3人以上でもよい）。
name date begin end
1 新井 2015/04/27 08:30 17:00
2 新井 2015/04/28 09:00 16:30
3 新井 2015/04/29 08:00 17:00
4 新井 2015/04/30 09:00 16:30
5 新井 2015/04/30 17:00 21:00
6 新井 2015/05/01 08:00 16:30
7 新井 2015/05/01 18:00 21:00
8 新井 2015/05/02 08:30 16:30
9 新井 2015/05/04 08:00 17:30
10 新井 2015/05/05 11:00 15:30
11 古田 2015/04/27 09:30 16:00
12 古田 2015/04/27 19:30 22:00
13 古田 2015/04/28 08:30 17:00
14 古田 2015/04/28 17:30 19:30
15 古田 2015/04/29 18:00 22:00
16 古田 2015/04/30 10:00 17:00
17 古田 2015/05/01 17:00 20:00
18 古田 2015/05/02 18:00 19:30
19 古田 2015/05/03 08:00 17:30
20 古田 2015/05/04 14:00 19:30
21 古田 2015/05/05 08:30 18:00
リストされている人全員が自由な時間が，連続する 3 時間以上となる，直近の日時を調べよ。ただし，6:00～22:00 の範囲であること。
func = function(s) {
# 6:00 ～ 22:00 までを 30 分を 1 コマとして何コマ目かを数字で表す
# 6:00～6:30 ==> 1, 6:30～7:00 ==> 2, ... , 21:30～22:00 ==> 32
hm = as.integer(unlist(strsplit(s, ":")))
hm[1] * 2 + (hm[2] == 30) - 11
}
inv = function(i) {
# 上記 func の逆関数コマ数を表す数字から，コマの開始時間を求める
sprintf("%02i:%02i", (i-1)%/%2+6, (i-1)%%2*30)
}
d = data.frame(stringsAsFactors = FALSE, # スケジュール表（文字列で）
name = c("新井", "新井", "新井", "新井", "新井", "新井", "新井",
"新井", "新井", "新井", "古田", "古田", "古田", "古田", "古田", "古田",
"古田", "古田", "古田", "古田", "古田"),
date = c("2015/04/27", "2015/04/28", "2015/04/29", "2015/04/30",
"2015/04/30", "2015/05/01", "2015/05/01", "2015/05/02",
"2015/05/04", "2015/05/05", "2015/04/27", "2015/04/27",
"2015/04/28", "2015/04/28", "2015/04/29", "2015/04/30",
"2015/05/01", "2015/05/02", "2015/05/03", "2015/05/04",
"2015/05/05"),
begin = c("08:30", "09:00", "08:00", "09:00", "17:00", "08:00",
"18:00", "08:30", "08:00", "11:00", "09:30", "19:30", "08:30",
"17:30", "18:00", "10:00", "17:00", "18:00", "08:00", "14:00",
"08:30"),
end = c("17:00", "16:30", "17:00", "16:30", "21:00", "16:30",
"21:00", "16:30", "17:30", "15:30", "16:00", "22:00", "17:00",
"19:30", "22:00", "17:00", "20:00", "19:30", "17:30", "19:30",
"18:00"))
d = d[order(d[, 2]), ] # 月日で並べ替え
n = nrow(d) # スケジュールの総数（行数）
t.begin = sapply(d[, 3], func) # 拘束開始時間をコマ数に変換
t.end = sapply(d[, 4], func) - 1 # 拘束終了時間をコマ数に変換
begin = 1 # 1 日分がどこからどこまでか
day = d[begin, 2]
for (i in 1:(n - 1)) {
if (d[i, 2] != d[i + 1, 2]) {
begin = c(begin, i + 1)
}
}
end = c(begin[-1] - 1, n)
for (i in seq_along(begin)) { # 日付ごとにチェック
time.table = integer(32) # コマ
for (j in begin[i]:end[i]) { # 拘束コマを 1 で埋める
time.table[t.begin[j]:t.end[j]] = 1
}
k = 1 # 連続する自由なコマ数を求める
while (k < 32) {
while (time.table[k] == 1) {
k = k + 1
if (k > 32) break
}
ks = k
if (k > 32) break
while (time.table[k] == 0) {
k = k + 1
if (k > 32) break
}
free.time = (k - ks) / 2 # コマ数を時間に直す（2 コマで 1 時間）
if (free.time >= 3) { # 自由時間が 3 時間以上なら結果出力
cat(sprintf("%s の %s ～ %s の %.1f 時間\n", d[j, 2], inv(ks), inv(k), free.time))
break
}
}
if (free.time >= 3) break
}
実行結果
2015/04/30 の 06:00 ～ 09:00 の 3.0 時間
(ai, bi), i = 1,2,...n において，
その部分集合 (am, bm), 1 ≦ m ≦ n を考える。
部分集合の全ての要素において j > k のときに aj > ak && bj > bk が成り立つ最大の部分集合の個数を求める。
(166, 71)
(178, 84)
(176, 94)
(174, 65)
(174, 65)
のときには，
(166, 71)
(174, 85)
(176, 94)
の 3 個の部分集合が求める解になる。
func = function(x) {
y = x[order(x[,1], x[,2]),]
n = nrow(y)
nrow.ans = 0
for (i in 1:(n-1)) {
ans = matrix(y[i,], nrow=1)
min1 = y[i, 1]
min2 = y[i, 2]
for (j in (i+1):n) {
if (y[j, 1] > min1 && y[j, 2] > min2) {
ans = rbind(ans, y[j,])
min1 = y[j, 1]
min2 = y[j, 2]
}
}
if (nrow(ans) > nrow.ans) {
nrow.ans = nrow(ans)
final.ans = ans
}
}
print(nrow.ans)
print(final.ans)
}
func(d1) # 3 個
func(d2) # 6 個
func(d3) # 37 個
d1, d2, d3 は，前もって以下のようにデータを入力しておく。
d1 = structure(c(166,178,176,174,174,71,84,94,85,65),.Dim = c(5,2),.Dimnames = list(NULL,c("V1","V2")))
d2 = structure(c(170,159,160,179,179,154,179,151,160,181,166,179,176,171,170,175,157,164,159,180,174,154,186,
161,168,156,183,178,158,158,74,72,65,92,46,61,77,61,49,49,62,48,50,65,69,65,60,55,58,56,54,54,82,47,62,50,67,
89,56,34),.Dim = c(30,2),.Dimnames = list(NULL,c("V1","V2")))
d3 = structure(c(182,180,164,171,170,187,161,163,171,180,158,183,163,173,151,166,167,171,175,154,155,177,154,
183,186,177,167,153,170,185,156,176,181,174,172,165,160,158,152,178,184,151,176,166,182,178,166,181,158,185,
174,183,177,164,178,162,158,152,155,153,156,172,155,160,173,162,168,185,158,166,153,171,186,155,160,188,163,
175,177,157,163,178,176,157,183,151,185,154,181,170,167,171,165,165,171,174,169,187,158,168,169,177,156,182,
162,166,163,189,188,169,151,165,158,188,188,170,187,156,157,168,157,176,180,154,158,158,186,172,153,169,162,
176,185,155,161,183,158,154,155,189,181,169,173,154,182,151,161,169,184,166,171,188,159,189,159,156,183,176,
175,164,171,156,184,184,177,156,150,180,150,177,179,180,185,174,179,179,161,174,159,167,164,151,150,183,182,
160,162,178,165,180,184,174,184,158,163,159,166,153,159,181,177,167,185,156,182,159,167,177,155,157,170,176,
188,180,167,158,154,176,171,171,167,154,173,173,178,185,171,171,173,182,165,189,165,163,185,176,161,176,187,
155,154,188,160,171,174,166,176,159,151,169,157,179,170,162,176,179,161,153,162,171,150,171,152,162,162,150,
185,151,183,151,170,158,188,177,185,166,172,159,168,189,152,172,184,163,153,168,150,178,186,177,160,166,170,
154,158,177,153,184,169,188,162,162,188,174,170,182,183,178,185,179,173,152,166,173,188,152,184,162,162,166,
185,173,161,166,169,163,168,183,165,157,165,183,175,166,154,172,167,154,151,188,167,189,186,171,172,188,167,
189,154,176,171,186,172,180,170,161,180,162,187,186,150,165,159,172,163,161,182,182,172,183,182,156,180,167,
159,174,173,171,175,177,160,157,186,185,181,160,178,162,173,179,172,173,188,153,174,175,181,183,164,182,163,
177,163,169,182,157,173,164,176,180,164,166,161,166,157,156,179,174,151,153,167,158,171,175,185,155,185,173,
176,156,159,172,154,187,169,186,158,177,186,174,164,168,186,184,159,155,175,189,150,155,189,163,167,172,185,
155,174,166,174,165,152,188,186,160,160,178,154,176,182,154,165,160,186,169,155,185,170,164,161,151,180,166,
176,170,179,179,150,181,167,167,173,170,169,155,155,171,155,153,163,166,159,158,178,188,166,189,173,179,181,
189,185,171,178,166,178,185,154,174,185,167,152,150,154,185,172,178,154,187,175,181,187,153,187,176,167,187,
177,151,153,175,181,184,159,183,152,180,168,164,166,178,176,154,187,166,167,151,185,175,176,166,181,189,175,
188,168,170,174,165,169,151,178,154,177,157,150,172,152,170,181,179,164,183,170,184,170,163,170,159,183,165,
151,160,169,183,154,171,177,175,173,163,171,175,186,174,177,167,180,179,154,150,188,175,158,184,187,183,161,
177,180,185,155,180,179,180,174,170,175,179,174,185,177,150,161,150,154,157,156,186,159,151,186,167,186,157,
160,186,159,179,178,165,168,174,174,184,166,177,155,154,179,172,187,178,170,184,154,169,154,161,182,166,181,
187,180,153,152,183,154,159,163,176,153,170,164,176,186,189,169,178,153,179,153,156,172,163,164,179,160,158,
182,186,181,177,155,160,154,171,150,174,171,159,158,181,156,164,171,170,162,181,159,164,156,175,179,187,151,
154,170,156,150,177,189,185,168,157,159,154,155,171,188,159,173,153,183,172,158,174,172,151,151,182,188,153,
157,167,178,152,168,154,187,163,160,161,171,150,183,174,188,186,185,162,162,153,175,184,152,176,176,150,162,
150,178,159,163,167,187,184,171,181,160,188,173,177,183,168,163,175,185,178,167,181,163,162,181,172,185,184,
177,182,150,178,154,185,185,161,166,176,166,155,184,174,187,153,158,170,168,160,150,163,153,150,182,152,163,
183,161,163,154,158,159,169,187,168,155,159,165,174,185,189,159,160,188,174,175,166,163,186,182,156,176,179,
168,169,156,173,171,161,161,169,154,181,156,163,170,167,164,172,185,185,181,153,186,182,185,180,183,151,182,
174,151,151,172,180,187,177,160,174,172,187,187,189,170,152,186,175,163,171,184,150,174,179,177,187,153,162,
151,161,174,189,150,181,164,154,158,153,166,155,165,187,184,168,179,161,185,185,183,162,177,163,164,180,166,
187,177,163,166,176,166,161,179,163,183,170,174,163,174,186,150,189,162,186,170,171,173,154,177,155,157,171,
161,165,175,150,164,171,176,174,162,177,173,180,162,156,188,171,186,155,168,189,177,165,160,183,188,167,152,
182,185,162,162,158,57,51,74,67,78,58,53,71,72,89,75,82,64,77,46,80,70,51,48,71,59,75,71,57,73,89,79,61,44,
72,67,71,62,56,72,62,53,75,56,50,70,36,50,74,87,63,51,59,62,54,71,67,51,66,59,64,35,55,57,48,70,39,61,59,44,
51,80,94,62,67,30,44,60,55,60,87,62,55,79,32,64,48,92,59,99,33,67,54,82,76,53,39,73,53,50,63,45,51,44,80,50,
45,54,97,37,59,80,101,85,86,43,52,55,101,78,83,84,33,72,80,39,91,66,68,70,65,55,79,58,79,37,55,104,43,67,72,
39,72,32,106,73,45,73,68,75,41,65,67,90,64,58,77,41,55,42,71,46,89,80,59,79,72,99,65,79,35,32,89,35,62,82,58,
55,43,71,44,75,90,61,66,78,55,34,50,66,53,52,77,63,78,70,77,78,37,61,46,61,35,50,43,48,37,60,48,59,70,61,71,
57,62,70,57,93,94,78,69,56,53,85,84,82,38,39,55,70,74,77,54,71,57,74,103,46,36,45,67,70,88,54,68,58,76,77,49,
73,37,42,41,44,76,32,50,53,56,85,64,39,59,68,61,68,55,62,47,60,58,92,59,72,62,54,51,85,69,56,74,72,46,81,57,
70,39,86,57,53,48,34,57,60,90,56,52,61,39,43,69,60,64,49,78,40,47,57,87,78,69,72,69,97,59,55,55,40,53,79,64,
100,71,56,51,64,55,47,47,42,60,82,93,80,46,68,73,90,44,52,70,76,40,51,62,71,55,91,58,75,57,67,99,48,57,59,95,
64,81,53,48,99,73,70,84,36,68,73,83,68,54,88,94,67,71,62,37,54,52,72,41,79,82,50,87,45,36,90,94,82,70,82,46,
55,81,56,70,72,40,55,88,59,88,40,84,37,73,39,53,91,57,63,41,53,44,53,63,39,63,71,62,64,87,41,35,84,49,70,46,
51,53,99,66,49,63,57,52,51,47,55,81,41,75,95,82,42,55,45,52,59,62,78,104,36,65,94,42,50,39,60,68,85,64,67,58,
35,84,51,43,41,73,35,50,75,45,78,77,46,83,65,81,48,41,70,46,48,50,87,46,60,93,47,73,77,78,43,57,63,57,67,46,
39,57,66,63,55,51,59,66,40,95,67,68,99,78,70,78,63,46,48,82,53,61,49,47,39,29,68,93,69,82,61,46,83,52,83,64,
80,71,39,90,56,42,63,92,81,48,61,86,41,58,61,68,60,60,93,48,102,46,83,45,94,55,65,41,48,52,90,71,48,66,67,59,
60,68,69,41,84,67,64,83,36,63,68,81,51,69,61,87,48,62,62,52,97,44,66,72,47,81,32,79,84,47,86,40,68,50,85,63,
84,42,96,90,38,38,94,51,36,93,83,61,65,76,52,102,34,93,79,79,61,60,62,78,55,94,90,30,41,64,62,49,41,63,73,47,
54,38,101,63,52,81,70,53,50,77,40,78,75,64,46,65,56,68,46,85,73,46,53,75,36,60,56,53,56,68,87,56,86,31,49,90,
59,72,40,66,46,51,68,79,92,98,60,84,35,78,52,39,61,39,53,47,73,64,96,75,60,77,50,72,56,59,52,89,75,75,52,95,
55,74,47,64,42,81,72,79,43,42,81,70,33,43,40,53,40,61,82,99,65,60,62,41,61,75,80,56,49,59,52,50,37,75,48,53,
42,44,108,54,59,67,68,67,62,42,100,43,65,44,41,67,69,88,55,56,101,58,52,42,44,58,62,51,57,43,59,63,68,49,39,
50,105,102,39,82,45,46,90,67,53,76,53,72,66,50,65,98,35,74,44,83,61,46,84,89,68,96,40,47,80,78,73,90,79,34,
51,43,65,50,73,71,50,38,60,47,35,56,66,47,66,63,45,54,49,73,40,72,81,43,60,53,37,70,98,62,65,52,97,48,57,65,
45,59,46,47,57,88,77,79,70,51,39,34,65,77,53,51,35,56,44,60,63,85,102,89,90,70,92,94,71,54,86,63,80,88,55,61,
82,95,103,84,76,61,46,66,73,62,59,69,64,49,73,86,96,64,63,57,84,74,57,38,48,75,71,91,30,53,82,36,34,52,78,32,
38,67,82,44,85,67,77,82,57,60,69,35,63,91,43,79,89,58,52,86,72,46,45,78,83,52,87,61,81,93,37,88,58,70,72,76,
40,43,51,53,53,65,72,46,93,57,37,63,74,47,56,84,74,79,69,65,52,39,97,47,86,84,46,76,55,49,53,65,58,53,86,55,
72,38),.Dim = c(1000,2),.Dimnames = list(NULL,c("V1","V2")))
2 つの時刻の間の時間数の筆算で，0 ～ 9 の全てが現れる計算は何通りあるか
例
16時27分
- 5時38分
==========
10時49分
0x 時とか 0y 分のような zero fill の 0 は除く
また，0 時 zz 分も除く
プログラム1
以下の AWK(GAWK) プログラムが 4 秒ほどで，一番速かった。
#!/usr/local/bin/gawk -f
BEGIN {
for (h1 = 1; h1 < 24; h1++) {
for (m1 = 0; m1 < 60; m1++) {
t1 = h1*60+m1
for (h2 = 1; h2 <= h1; h2++) {
for (m2 = 0; m2 < 60; m2++) {
t2 = h2*60+m2
if (t2 > t1) break
diff = t1-t2
m3 = diff%60
h3 = int(diff/60)
tbl(h1)
tbl(h2)
if (h3 > 0) tbl(h3)
tbl(m1)
tbl(m2)
tbl(m3)
if (length(digits) == 10) {
# print h1, m1, h2, m2, h3, m3
count++
}
delete digits
}
}
}
}
print "ans =", count
}
function tbl(d, d1) {
d1 = int(d/10)
if (d1 > 0) digits[d1]
digits[d%10]
}
文字列を使って以下のようにしたものは 7 秒強かかる。
#!/usr/local/bin/gawk -f
BEGIN {
for (h1 = 1; h1 < 24; h1++) {
for (m1 = 0; m1 < 60; m1++) {
t1 = h1*60+m1
for (h2 = 1; h2 <= h1; h2++) {
for (m2 = 0; m2 < 60; m2++) {
t2 = h2*60+m2
if (t2 > t1) break
diff = t1-t2
m3 = diff%60
h3 = int(diff/60)
hm = m1 m2 m3 h1 h2
if (h3 > 0) hm = hm h3
split(hm, a, "")
for (i in a) digits[a[i]]
if (length(digits) == 10) {
# print h1, m1, h2, m2, h3, m3
count++
}
delete digits
}
}
}
}
print "ans =", count
}
最初のプログラムを R で描き直すと，25 秒もかかる。ほぼ 6 倍ということだ。
tbl = function(d, digits) {
d1 = d %/% 10
if (d1 > 0) digits[d1+1] = TRUE
digits[d %% 10+1] = TRUE
return(digits)
}
count = 0
for (h1 in 1:23) {
for (m1 in 0:59) {
t1 = h1*60+m1
for (h2 in 1:h1) {
for (m2 in 0:59) {
t2 = h2*60+m2
if (t1 < t2) break
diff = t1-t2
m3 = diff %% 60
h3 = diff %/% 60
digits = logical(10)
digits = tbl(h1, digits)
digits = tbl(h2, digits)
if (h3 > 0) digits = tbl(h3, digits)
digits = tbl(m1, digits)
digits = tbl(m2, digits)
digits = tbl(m3, digits)
if (sum(digits) == 10) {
# cat(h1, m1, h2, m2, h3, m3, "\n")
count = count+1
}
}
}
}
}
cat("ans =", count, "\n")
R で定石のベクトル演算を使って書いてみたら，260 秒もかかるとんでもないプログラムになった。最短 AWK プログラムのほぼ 65 倍ということだ。
a = expand.grid(1:23, 0:59, 1:23, 0:59)
t1 = a[,1]*60+a[,2]
t2 = a[,3]*60+a[,4]
diff = t1 - t2
flag = diff > 0
a$diff1 = diff %/% 60
a$diff2 = diff %% 60
a = a[flag,]
func = function(d) {
d2 = d%%10
d1 = d%/%10
d1[d1 == 0] = NA
return(length(table(c(d1, d2))))
}
count = 0
apply(a, 1, function(x) {
if (x[5] == 0) x = x[-5]
if (func(x) == 10) {
# cat(h1, m1, h2, m2, h3, m3, "\n")
count+1 ->> count
}
})
cat("ans =", count, "\n")
解答は，以下の 48 通りということである。
16時27分 - 5時38分 = 10時49分
16時27分 - 5時39分 = 10時48分
16時27分 - 5時48分 = 10時39分
16時27分 - 5時49分 = 10時38分
16時27分 - 10時38分 = 5時49分
16時27分 - 10時39分 = 5時48分
16時27分 - 10時48分 = 5時39分
16時27分 - 10時49分 = 5時38分
18時25分 - 7時36分 = 10時49分
18時25分 - 7時39分 = 10時46分
18時25分 - 7時46分 = 10時39分
18時25分 - 7時49分 = 10時36分
18時25分 - 10時36分 = 7時49分
18時25分 - 10時39分 = 7時46分
18時25分 - 10時46分 = 7時39分
18時25分 - 10時49分 = 7時36分
20時47分 - 3時48分 = 16時59分
20時47分 - 3時49分 = 16時58分
20時47分 - 3時58分 = 16時49分
20時47分 - 3時59分 = 16時48分
20時47分 - 6時48分 = 13時59分
20時47分 - 6時49分 = 13時58分
20時47分 - 6時58分 = 13時49分
20時47分 - 6時59分 = 13時48分
20時47分 - 13時48分 = 6時59分
20時47分 - 13時49分 = 6時58分
20時47分 - 13時58分 = 6時49分
20時47分 - 13時59分 = 6時48分
20時47分 - 16時48分 = 3時59分
20時47分 - 16時49分 = 3時58分
20時47分 - 16時58分 = 3時49分
20時47分 - 16時59分 = 3時48分
20時57分 - 4時18分 = 16時39分
20時57分 - 4時19分 = 16時38分
20時57分 - 4時38分 = 16時19分
20時57分 - 4時39分 = 16時18分
20時57分 - 6時18分 = 14時39分
20時57分 - 6時19分 = 14時38分
20時57分 - 6時38分 = 14時19分
20時57分 - 6時39分 = 14時18分
20時57分 - 14時18分 = 6時39分
20時57分 - 14時19分 = 6時38分
20時57分 - 14時38分 = 6時19分
20時57分 - 14時39分 = 6時18分
20時57分 - 16時18分 = 4時39分
20時57分 - 16時19分 = 4時38分
20時57分 - 16時38分 = 4時19分
20時57分 - 16時39分 = 4時18分
R だと，
> ×
>
> s = 0
>
> for (n in 0:9)
> {
> s = s+n
> }
>
> print(s)
>
>
>
> ◯
>
> s = 0:9
>
> Reduce("+", s)
sum(s) で十分じゃないのかな？
sum の中で何がどう扱われていようと，見かけは関数に過ぎないわけで。
主成分分析だって prcomp 関数一つでできることになるんだしね。その中では「フローチャート」によるアルゴリズムの実体化が行われているし，計算もループによる収束計算を含んだりしているわけで。
フィボナッチ数列の各項には，行列を用いて以下のような関係がある。
この行列を活用して，数同士の足し算や掛け算の回数を 200 回以内に抑えてフィボナッチ数列の第 1000 項の下 5 桁を求めよ。
解答としては「フィボナッチ数列の第 1000 項の下 5 桁」と「計算中の数同士の足し算や掛け算の回数」を答えよ。
ただし、数を下 5 桁にする計算は足し算や掛け算として含めない。
行列を t とすると，t の n 乗 の 2 行目の 1 列目の要素がフィボナッチ数列の n 項目になる。
> a = t = matrix(c(1,1,1,0), 2)
> t
[,1] [,2]
[1,] 1 1
[2,] 1 0
> for (i in 2:10) a = a %*% t
> a
[,1] [,2]
[1,] 89 55
[2,] 55 34
> library(gmp)
> fibnum(10)
Big Integer ('bigz') :
[1] 55
1000 項目は t の 1000 乗を求めればよいが，t の 512, 256, 128, 64, 32, 8 乗 の積を求めるのが効率がよい。
なお，フィボナッチ数列の末尾 5 桁を求めればよいので，行列の積を求めるごとに，行列の要素を 100000 で割った余りを演算結果として残せばよい。
t = matrix(c(1,1,1,0), 2)
t2 = t %*% t %% 100000
t4 = t2 %*% t2 %% 100000
t8 = t4 %*% t4 %% 100000
t16 = t8 %*% t8 %% 100000
t32 = t16 %*% t16 %% 100000
t64 = t32 %*% t32 %% 100000
t128 = t64 %*% t64 %% 100000
t256 = t128 %*% t128 %% 100000
t512 = t256 %*% t256 %% 100000
a = t512 %*% t256 %% 100000
a = a %*% t128 %% 100000
a = a %*% t64 %% 100000
a = a %*% t32 %% 100000
a = a %*% t8 %% 100000
> a
[,1] [,2]
[1,] 3501 28875
[2,] 28875 74626
1000 項目の末尾 5 桁は a[2,1] で 28875 である。
単純に 3,4,5,..., 1000 項を求めるプログラムは以下のようになる。足し算は 998 回だ。
b = integer(1000)
b[1] = b[2] = 1
for (i in 3:1000) b[i] = (b[i-2]+b[i-1]) %% 100000
b[1000]
> b[1000]
[1] 28875
同じ答えが得られることが確かめられた。
2×2 行列の積は，結果の各要素ごとに 2 回の掛け算と，1 回の足し算が必要である。上記プログラムは 14 回の行列積を行う。
よって，全ての計算において 14 * (2 * 4) = 112 回の掛け算と，14 * (1 * 4) = 56 回の足し算が行われる。
998 回の足し算と
112 回の掛け算と 56 回の足し算
どちらがコストが高いか。
> fib1 = function() {
+ t = matrix(c(1,1,1,0), 2)
+ t2 = t %*% t %% 100000
+ t4 = t2 %*% t2 %% 100000
+ t8 = t4 %*% t4 %% 100000
+ t16 = t8 %*% t8 %% 100000
+ t32 = t16 %*% t16 %% 100000
+ t64 = t32 %*% t32 %% 100000
+ t128 = t64 %*% t64 %% 100000
+ t256 = t128 %*% t128 %% 100000
+ t512 = t256 %*% t256 %% 100000
+ a = t512 %*% t256 %% 100000
+ a = a %*% t128 %% 100000
+ a = a %*% t64 %% 100000
+ a = a %*% t32 %% 100000
+ a = a %*% t8 %% 100000
+ }
>
> fib2 = function() {
+ b = integer(1000)
+ b[1] = b[2] = 1
+ for (i in 3:1000) b[i] = (b[i-2]+b[i-1]) %% 100000
+ b[1000]
+ }
>
> library(rbenchmark)
> benchmark(fib1, fib2, replications=100000)
test replications elapsed relative user.self sys.self user.child sys.child
1 fib1 100000 0.417 1.185 0.364 0.040 0 0
2 fib2 100000 0.352 1.000 0.348 0.007 0 0
ほとんど差はない，わざわざ複雑な計算方法を採る意味がないということになった。（最初から分かっている？）
PVアクセスランキング にほんブログ村
