覆面算

READ + WRITE + TALK = SKILL　を満たす解が何通りあるか。

以下のプログラムは，2 番目のものより 10 倍ほど「遅い」。

library(e1071)
d = permutations(10)-1
w = c(1000, 100, 10, 1)
w2 = c(w*10, 1)
ans = apply(d, 1, function(x) all(x[c(1,5,7,10)] != 0) && sum(x[1:4]*w)+sum(x[c(5,1,6,7,2)]*w2)+sum(x[c(7,3,8,9)]*w) == sum(x[c(10,9,6,8,8)]*w2))
sum(ans)
d = d[ans,]
apply(d, 1, function(x) cat(sprintf("%i + %i + %i = %i\n", sum(x[1:4]*w), sum(x[c(5,1,6,7,2)]*w2), sum(x[c(7,3,8,9)]*w), sum(x[c(10,9,6,8,8)]*w2))))

なるべくベクトル演算を使うと速くなる。

library(e1071)
d = t(permutations(10)-1)
w = c(1000, 100, 10, 1)
w2 = c(w*10, 1)
read  = colSums(d[1:4,]*w)
write = colSums(d[c(5,1,6,7,2),]*w2)
talk  = colSums(d[c(7,3,8,9),]*w)
skill = colSums(d[c(10,9,6,8,8),]*w2)
ans = d[1,]*d[5,]*d[7,]*d[10,] != 0 & read+write+talk == skill
sum(ans)
data.frame(read=read[ans], write=write[ans], talk=talk[ans], skill=skill[ans])

   read write talk skill
1  5270 85132 3764 94166
2  2543 72065 6491 81099
3  5180 65921 2843 73944
4  1632 41976 7380 50988
5  5094 75310 1962 82366
6  5096 35710 1982 42788
7  7092 37510 1986 46588
8  7092 47310 1986 56388
9  4905 24689 8017 37611
10 9728 19467 6205 35400

 

添え字の検索

プログラムの冒頭に示すような，1 文字を要素とする 20×20 の行列がある。

"STAY HUNGRY, STAY FOOLISH" の各文字が，行列のどこにあるか，添え字を求めよ。

ただし，一度検索された要素（添え字）は二度目以降は検索対象にならない。

box <- matrix(c(
"Z","T","A","V","O","J","Z","G","W","L","M","A","W","B","Y","I","M","D","K","Y",
"P","G","P","N","C","Q","E","R","I","Q","I","T","O","K","H","F","T","O","Z","Z",
"D","T","S","E","L","F","N","O","L","H","Z","K","L","N","F","D","J","C","G","A",
"A","E","U","J","U","I","N","Z","D","T","D","U","W","J","O","N","C","S","L","H",
"Z","W","E","F","G","Z","B","R","I","O","T","U","N","B","O","G","N","U","Z","Y",
"U","Y","C","O","C","P","Z","G","W","C","O","D","J","X","J","L","F","V","J","H",
"D","B","K","Q","I","J","M","Q","X","P","U","P","Z","G","C","G","S","M","I","P",
"F","E","N","Q","R","Z","T","B","U","J","T","M","E","B","U","T","N","F","E","Y",
"R","E"," ","Y","D","G","B","A","C","F","B"," ","A","B","W","X","U","J","A","Z",
"M","D","O","W","W","X","S","C","O","E","J","W","N","L","U","F","E","X","H","X",
"E","F","O","F","W","P","V","R","U","I","G","R","E","Z","W","O","U","G"," ","E",
"W","I","M","V","D","Q","V","S","I","Y","Y","B","R","Z","O","Y","S","U","W","Y",
"D",",","A","G","U","F","L","Z","L","N","X","O","E","Q","J","X","E","G","H","E",
"O","N","E","M","V","J","Z","N","C","Q","P","Z","K","W","F","M","V","K","Y","M",
"M","B"," ","H","R","M","L","S","T","I","A","O","O","O","O",",","E","H","M","E",
"I","U","O","Q","B","J","Z","N","Q","N","Q","H","U","D","Y","J","A","B","I","P",
"X","K","K","I","S","H","Z","H","D","L","Z","O","J","N","P","C","W","O","O","R",
"L","Y","Y","M"," ","B","N","T","K","X","C"," ","S","B","K","H","R","W","M","D",
"V","N","U","L","L","Z","M","I","O","G","S","M","P","U","E","P","E","O","A","C",
"R","E","T","L","F","S","D","T","Q"," ","U","P","Y","K","A","F","P","G","W","O"),
ncol=20, byrow=TRUE)
types = unlist(strsplit("STAY HUNGRY, STAY FOOLISH", ""))
for (i in seq_along(types)) {
  ans = which(types[i] == box, arr.ind=TRUE)[1, , drop=FALSE]
  box[ans] = ""
  print(c(ans))
}

最適な組合せ

16 種類のオーダーと，それにかかる日数と売上代金の表がある。同じオーダーを 2 つ以上使うことはできない。同時に 2 つ以上のオーダーは処理できない。代金が最大になるような複数のオーダーの組合せを求めよ。ただし，所要日数の合計は 100 日以内でなければならない。

            所用日数  売上代金
オーダー 1         9      130
オーダー 2        12      150
オーダー 3        20      190
オーダー 4        23      190
オーダー 5        27      230
オーダー 6        33      290
オーダー 7        31      330
オーダー 8         9       70
オーダー 9        30      330
オーダー 10        9      110
オーダー 11        6       90
オーダー 12       34      310
オーダー 13       34      330
オーダー 14       22      190
オーダー 15       25      230
オーダー 16       13      170

fee = rev(c(130, 150, 190, 190, 230, 290, 330, 70, 330, 110, 90, 310, 330, 190, 230, 170))
day = rev(c(9, 12, 20, 23, 27, 33, 31, 9, 30, 9, 6, 34, 34, 22, 25, 13))
library(e1071)
bc = bincombinations(length(fee))
ans = apply(bc, 1, function(x) { x = as.logical(x); ifelse(sum(day[x]) <= 100, sum(fee[x]), 0) })
( max.fee = max(ans) )
apply(bc[ans == max.fee, ], 1, function(x) paste(paste("Order", which(rev(x) == 1)), collapse=", "))

最短距離のルート選択

「最もエネルギー効率のよいルート選択」と同じであるが，こちらは，合計距離が最短になるルート選択。エネルギー行列と違って，距離行列は対称行列である。

ノードが少なければ，しらみつぶし探索でも問題ないが，ノードが多くなったら，別の解法を採らないとね...

dist = matrix(c(0,1,3,4,5, 1,0,2,3,4, 3,2,0,1,2, 4,3,1,0,1, 5,4,2,1,0), 5)
library(e1071)
route = permutations(5)
d = apply(route, 1, function(r) sum(dist[cbind(r, c(r[-1], r[1]))]))
( min.dist = min(d) ) # 最短距離
route[d == min.dist,] # 径路（逆順も含む）

最もエネルギー効率のよいルート選択

6 つのポイントがあり，あるポイントから別のポイントへ移動するために必要なエネルギーが energy.matrix で記述されている。1 番目のポイントから 2 番目のポイントに移動するためには，energy.matrix[1, 2] = 7 必要というようなこと。

1 番目のポイントからスタートして，全てのポイントを回った後，最後に 1 番目のポイントに帰る。どのようなルートをとれば，最小エネルギーで回れるか。

解答例

n = 6
energy.matrix <- matrix(c(
0,7,12,8,11,7,
3,0,10,7,13,2,
4,8,0,9,12,3,
6,6,9,0,10,7,
7,7,11,10,0,5,
9,7,8,9,10,0),
ncol=n, byrow=TRUE)
library(e1071)
route = permutations(n)
route = route[route[,1] == 1,] # 1 番目からスタートするルートのみ
total.energy = apply(route, 1, function(r) sum(energy.matrix[cbind(r, c(r[-1], 1))]))
where = which.min(total.energy)
route[where,]
total.energy[where]

実行例

> route[where,] # ルート
[1] 1 4 5 2 6 3

> total.energy[where] # そのときのエネルギー
[1] 39
    
> sum(total.energy == 39) # 最小値をとるルートはいくつあるか確認
[1] 1

> table(total.energy) # 必要エネルギーの分布
total.energy
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 59
 1  1  3  5  8 10  8 13 15 11 12 10  9  3  5  2  2  1  1

Color Space の変換

特に難しいものではないと思う。

プログラム検証に例が100個も必要か？

CS2BS = function(CS) {
  mx = matrix(c(0,1,2,3, 1,0,3,2, 2,3,0,1, 3,2,1,0), 4)
  name = c("A","C","G","T")
  dimnames(mx) = list(name, name)
  BS = unlist(strsplit(CS, ""))
  former = which(BS[1] == name)
  for (i in 2:nchar(CS)) {
    BS[i] = names(which(BS[i] == mx[former,]))
    former = which(BS[i] == name)
  }
  paste(BS, collapse="")
}

sample01 = "C103023011232302323231301122233320313322132310333233003233231023023000231223"
answer01 = "CAATTCGGTGATCGGATCGATGCCACTCTATAGGCATAGACGATGGCGCTATTTAGCGATGGATTCGGGGATGAGC"
CS2BS(sample01) == answer01

sample02 = "G013232302011113133323231113100032313121302131301303000310323130301323123032"
answer02 = "GGTAGCTAAGGTGTGCATATCGATGTGCAAAATCGTACTGCCTGCATTGCCGGGGCAATCGTAATTGCTACTAATC"
CS2BS(sample02) == answer02

5×5ラテン方陣アナグラム

数限られているから数回の試行錯誤でもできるとは思うが，ユニーク解であるかどうかは，しらみつぶしでやるしかないかな

word.list = c("steal", "stela", "telas", "teals", "elast", "least", "laste", "astel")
w.l = sapply(word.list, function(w) unlist(strsplit(w, "")))
ans = combn(length(word.list), nchar(word.list[1]), function(o) apply(w.l[,o], 1, paste, collapse=""))
sapply(ans[,apply(apply(ans, 2, "%in%", word.list), 2, all)], function(w) unlist(strsplit(w, "")))

combn では，全てを尽くせないので，e1071 の permutations を使って以下のようにやってみる。

word.list = c("steal", "stela", "telas", "teals", "elast", "least", "laste", "astel")
w.l = sapply(word.list, function(w) unlist(strsplit(w, "")))
suf = unique(permutations(length(word.list))[, 1:nrow(w.l)])
ans = apply(suf, 1, function(o) apply(w.l[,o], 1, paste, collapse=""))
ans = sapply(ans[,apply(apply(ans, 2, "%in%", word.list), 2, all)], function(w) unlist(strsplit(w, "")))
dim(ans) = c(5, 5, ncol(ans)/5)
ans

換字式暗号

エドガー・アラン・ポーの黄金虫より簡単だ（情報が多いからね）

R で書くともっと簡単。数回の試行錯誤で解読可能。

chartr(
"ryqzvngfluohticxdbamwpejks",
"etonasirhdlucfmpgjwyvkbzxq",
tolower(readLines("text.txt")))

整数 a から 整数 b までに数字 n は何個あるか

library(gmp)

count2 = function(s, n) {
  if (nchar(s) == 1) {
    return(as.bigz(as.integer(s) >= n))
  } else {
    left = as.integer(substring(s, 1, 1))
    s = substring(s, 2, nchar(s))
    result = left*nchar(s)*as.bigz(10)^(nchar(s)-1) -
             (n == 0)*as.bigz(10)^(nchar(s))
    if (left == n) {
      result = result + as.bigz(sub("^0+([0-9]+)", "\\1", s))+1
    } else if (left > n) {
      result = result + as.bigz(10)^nchar(s)
    }
    return(result + Recall(s, n))
  }
}

# 整数 a から 整数 b までに数字 n は何個あるか
# a, b は文字列，n は数値で与える
count = function(a, b, n) {
  n.b = count2(b, n)
  n.a = if (as.bigz(a) == 0) as.bigz(0) else count2(as.character(as.bigz(a)-1), n)
  n.b-n.a
}
 
> print(count("14", "72", 1), initLine=FALSE) # ans = 12 , ans0 = 2 FALSE
[1] 12
> print(count("193", "201", 9), initLine=FALSE) # ans = 8 , ans0 = 18 FALSE
[1] 8
> print(count("144", "218", 9), initLine=FALSE) # ans = 17 , ans0 = 27 FALSE
[1] 17
> print(count("558", "960", 9), initLine=FALSE) # ans = 142 , ans0 = 141 FALSE
[1] 142
> print(count("6211", "19382", 9), initLine=FALSE) # ans = 5310 , ans0 = 5300 FALSE
[1] 5310
> print(count("10416", "18923", 5), initLine=FALSE) # ans = 3600 , ans0 = 3601 FALSE
[1] 3600
> print(count("288292", "935605", 0), initLine=FALSE) # ans = 329167 , ans0 = 329168 FALSE
[1] 329167
> print(count("633636", "635889", 2), initLine=FALSE) # ans = 645 , ans0 = 655 FALSE
[1] 645
> print(count("151642695", "184774424", 0), initLine=FALSE) # ans = 22862553 , ans0 = 22862552 FALSE
[1] 22862553
> print(count("52586824", "348090208", 0), initLine=FALSE) # ans = 236341495 , ans0 = 236341485 FALSE
[1] 236341495
> print(count("135279539884", "865804912485", 1), initLine=FALSE) # ans = 865138642422 , ans0 = 865138642421 FALSE
[1] 865138642422
> print(count("1761231736987", "22345608831150", 0), initLine=FALSE) # ans = 26743902079488 , ans0 = 26743902079487 FALSE
[1] 26743902079488
> print(count("6535079026811", "20167019383277", 9), initLine=FALSE) # ans = 18391758529453 , ans0 = 18391758529454 FALSE
[1] 18391758529453

CodeIQ 7 を数える（2）

「CodeIQ 7 を数える」を任意の数字の個数を数えることができるように，再帰関数として定義した。

count2 = function(s, n) {
  if (nchar(s) == 1) {
    return((as.integer(s) >= n)+0)
  } else {
    left = as.integer(substring(s, 1, 1))
    s = substring(s, 2, nchar(s))
    count = left*nchar(s)*10^(nchar(s)-1) - (n == 0)*10^(nchar(s))
    if (left == n) {
      count = count + as.integer(s)+1
    } else if (left > n) {
      count = count + 10^nchar(s)
    }
    return(count + Recall(s, n))
  }
}

count = function(s, n) {
  sum(unlist(strsplit(as.character(0:as.integer(s)), "")) == n)
}

> count2("649", 1)
[1] 235
> count("649", 1)
[1] 235

> count2("3649", 0)
[1] 1025
> count("3649", 0)
[1] 1025

これをベースにして他の問題を解こう。

円周上に並んでいる個体を 1 個ずつ取り除いて...

1～2*n 番目までの 2*n 個の個体が円周上に並んでいる。2 番目から初めて 1 つ置きに取り去る。最後に取り除かれるのが最初に並んだときの n 番目の個体であるような n を小さい方から 5 つ求めよ。

k = 0
n = 5 # 21, 85, 341, 1365, 5461, 21845       i.e. a[1]=5, a[i]=4*a[i-1]+1, i > 1
repeat {
  n = n+1
  seat = logical(2*n)
  seat[n] = TRUE
  i = 2
  ok = TRUE
  for (m in (2*n-1):1) {
    if (seat[i]) {
      ok = FALSE
      break
    }
    seat
    i = i+1
    if (i > m) i = i-m
  }
  if (ok && seat) {
    print(n)
    k = k+1
    if (k == 5) break
  }
}

このプログラムで 4 つ目までは比較的短時間で求まる。最初の n=5 を含めて，よく眺めると，規則性があることが分かる。で，5 つ目は 5461 であると当たりを付けて計算してみると確かに n = 5461 は解である。21845 は 6 つ目...

証明ができればよいのだけど。

漸化式が知られている

> 12人全員が他人のパソコンを持って帰ってしまうような組み合わせが何通りあるか

Wikipedia 曰く，

完全順列（かんぜんじゅんれつ、英 Derangement）、もしくは攪乱順列（かくらんじゅんれつ）とは、整数 1, 2, 3, …, n を要素とする順列において、i 番目 (i ≤ n) が i でない順列である。

完全順列の総数はモンモール数という。

モンモール数 a_n を与える漸化式は a_n = (n-1)(a_n-1+a_n-2), n ≧ 3

a₁ = 0, a₂ = 1

> an = function(n) {
+     if (n == 1) {
+         return(0)
+     } else if (n == 2) {
+         return(1)
+     } else {
+         return((n-1)*(Recall(n-1)+Recall(n-2)))
+     }
+ }

> an(3)
[1] 2

> an(4)
[1] 9

> an(12)
[1] 176214841

変わった数の二乗

9 桁の整数を二乗して，結果の下 9 桁が最初の数と同じになる。そういう数を見つけなさいと...

素直にやっていくとどうしようもない。どうせ，各桁の取り得る値は限られているのだから，下の桁から可能性を絞っていくという戦略で，ループを段階的に深くしたプログラムにしていく。問題は 9 桁の数の二乗は倍精度実数でも整数として表現しきれない場合があるということ。最初から gmp を使うのも手ではあるが，そこはぐっと我慢して最後の桁のチェックは別法にて行う。

for (i in c(0, 1, 5, 6)) {
  for (h in c(0, 2, 7)) {
    for (g in c(3, 6)) {
      for (f in c(0, 9)) {
        for (e in c(0, 9)) {
          for (d in c(1, 8)) {
            for (c in c(2, 7)) {
              for (b in c(1, 8)) {
                for (a in 0:9) {
#                  x = ((((((((a*10+b)*10+c)*10+d)*10+e)*10+f)*10+g)*10+h)*10+i)^2
#                  s = as.integer(rev(unlist(strsplit(as.character(x), ""))))
                  s = integer(10)
                  s[1] = i*i
                  s[2] = 2*h*i
                  s[3] = 2*g*i+h*h
                  s[4] = 2*(f*i+g*h)
                  s[5] = 2*(e*i+f*h)+g*g
                  s[6] = 2*(d*i+e*h+f*g)
                  s[7] = 2*(c*i+d*h+e*g)+f*f
                  s[8] = 2*(b*i+c*h+d*g+e*f)
                  s[9] = 2*(a*i+b*h+c*g+d*f)+e*e
                  for (k in 1:9) {
                    if (s[k] > 9) {
                      s[k+1] = s[k+1] + s[k] %/% 10
                      s[k] = s[k] %% 10
                    }
                  }
                  if (length(s) >= 9 && s[1] == i && s[2] == h && s[3] == g && s[4] == f && s[5] == e && s[6] == d && s[7] == c && s[8] == b && s[9] == a) {
                    answer = paste(a, b, c, d, e, f, g, h, i, sep="")
                    print(answer, initLine=FALSE)
                    print(as.bigz(answer)^2)
                  }
                }
              }
            }
          }
        }
      }
    }
  }
}

データの平均値（合計）を求めるとき...

データの性質を見極めておかないと誤った統計数値を得てしまうよという教訓を垂れる訳だろうが...

図に示すような 3 種類のデータがある。各データは 2 列で，1 列目と 2 列目の合計が等しいかどうか判断せよと。

1 番目のデータ（左上が 1 列目，左下が 2 列目）のヒストグラム
　かなり大きな正の値のデータのみ

2 番目のデータ（中上が 1 列目，中下が 2 列目）のヒストグラム
　正・負のデータからなり，絶対値がほどほど大きいデータを少し含む

3 番目のデータ（右上が 1 列目，右下が 2 列目）のヒストグラム
　正・負のデータからなり，外れ値ともいえるような絶対値がめちゃ大きいデータをほんのわずか含む

無邪気に，mean を使ったりすると，

> cat(sum(q1[,1]), sum(q1[,2]))
9007199254740992 9007199254740992

> cat(sum(q2[,1]), sum(q2[,2]))
0 -4.54747350886464e-13

> cat(sum(q3[,1]), sum(q3[,2]))
-71606554540133104 8756878575033.26

のようになる。1 番目のデータの各列の和は等しく，2 番目のデータでは誤差範囲で同じ，3 番目のデータではすごく違いがあるように見えるかもしれないが元のデータのオーダーが 30 桁ほどなんだから，誤差範囲であろうと判断しなければならない。

そもそも，2 番目とか 3 番目のようなデータは実際にはない（オーダーが大きいのは，単位をかえればよいし，普通そのように表されるだろう）ので，そんなに神経質にならなくてもよいのだけど。

それでも，次の様なアルゴリズムで和を求めるというのは，予防的措置として推奨されるかもしれない。

a. 足し算と引き算を混ぜない。正のデータと負のデータを別々に和を求めて，最後に 2 つの和を合計する。

b. 小さな数と大きな数の足し算をしない。小さなものから順に足していく（より慎重にするなら，2 つのデータの和を求めたら，残りのデータと一緒にしてその中の小さな 2 つのデータの和を求めるということを繰り返す）。

以下のようなプログラムを書いてみる。

sum2 = function(d) {
    sum1 = function(x) {
        x1 = x[x > 0]
        x2 = x[x < 0]
        cat(sum(x1), sum(abs(x2)), sum(x1)-sum(abs(x2)), "\n")
        ifelse(length(x1) > 0, tail(cumsum(sort(x1)), 1), 0) -
        ifelse(length(x2) > 0, tail(cumsum(sort(abs(x2))), 1), 0)
    }
    sum1(d[,1]) == sum1(d[,2])
}

これによる結果は，以下のようになる。

> sum2(q1)
9007199254740992 0 9007199254740992
9007199254740992 0 9007199254740992
[1] TRUE

> sum2(q2)
196452.384765625 196452.384765625 0
16530519.4 16530519.4 0
[1] TRUE

> sum2(q3)
2.76019691920891e+33 2.76019691920891e+33 0
1.82694806300339e+29 1.82694806300339e+29 35184372088832
[1] FALSE

3 番目のデータにおいてはまだ誤差範囲であるが同じではないという結果である。

えい，まだるっこしい。

いっそのこと，長精度実数計算をしてしまえ（ただし，上に書いた a, b は考慮する）。

library(Rmpfr)
sum9 = function(d) {
    sum8 = function(x) {
        prec = 200
        x1 = mpfr(x[x > 0], prec)
        x2 = mpfr(x[x < 0], prec)
        s1 = if (length(x1) > 0) tail(cumsum(sort(x1)), 1) else 0
        s2 = if (length(x2) > 0) tail(cumsum(sort(abs(x2))), 1) else 0
        print(c(s1, s2, s1-s2))
        s1-s2
    }
    sum8(d[,1]) == sum8(d[,2])
}

これによれば，全てのデータセットで，2 列の合計値は等しくないということがはっきりする。

そして，その差は誤差範囲であるということである。

> sum9(q1)
3 'mpfr' numbers of precision  128 .. 200  bits
[1] 9007199254740993                0 9007199254740993
3 'mpfr' numbers of precision  128 .. 200  bits
[1] 9007199254740992                0 9007199254740992
[1] FALSE

> sum9(q2)
3 'mpfr' numbers of precision  200   bits
[1] 196452.384765625 196452.384765625                0
3 'mpfr' numbers of precision  200   bits
[1]         16530519.40000000000000479616346638067625463008880615234375
[2]        16530519.399999999999994582111639829236082732677459716796875
[3] 1.0214051826551440171897411346435546875000000000000000000000000e-14
[1] FALSE

> sum9(q3)
3 'mpfr' numbers of precision  200   bits
[1]  2760196919208909839959496867497081.8282129668101326662451169565
[2]  2760196919208909911700429409141880.8282129668101326104341159326
[3] -71740932541644798.999999999999999944188998976167179947377656232
3 'mpfr' numbers of precision  200   bits
[1] 182694806300338729143596181240.44380035052800129831732213181835
[2] 182694806300338720381823051512.44380035052800286368844502856631
[3] 8761773129727.9999999999999984346288771032520431883572850254938
[1] FALSE

 

PI には何でもありということで...

自分で PI の小数点以下の数字を求めて，その中に特有の数字列がどこにあるかを探すプログラムを書く。PI の数字列は自分で計算することを条件とする。

たとえば 9 が 6 個連続して出現するのはどこか？
そのほかに，どのような「珍しい」数字列がどこにあるかを探索するプログラムを書く。