RでProject Euler (Problem 51~60)

Problem 41~50

Problem 51 「素数の桁置換」

$\ast 3$の第1桁を置き換えることで,$13, 23, 43, 53, 73, 83$という6つの素数が得られる.
$56\ast\ast3$の第3桁と第4桁を同じ数で置き換えることを考えよう. この5桁の数は7つの素数をもつ最初の例である: $56003, 56113, 56333, 56443, 56663, 56773, 56993$. よって, この族の最初の数である$56003は$, このような性質を持つ最小の素数である.
桁を同じ数で置き換えることで8つの素数が得られる最小の素数を求めよ. (注:連続した桁でなくても良い)

・ 一の位が$\ast$だと、偶数や5に置き換えたときに素数になることはない
・ $\ast$が0か1か2のときに素数になる(例えば3が最小だと3~9で7つの素数しか得られない)
・ $\ast$の総数は、3の倍数個ある($\ast$の数が例えば1個や2個のとき、3の倍数を考えると最大でも7つの素数しか得られない)
以上の条件で絞りつつ探索する。

上限が見えないため、上限を10倍ずつして求める数を探す。

getprimes = function(limit){
  p = rep(TRUE, limit)
  p[1] = FALSE
  for(i in c(2, seq(3, sqrt(limit), 2)))
    if(p[i])
      p[seq(i*2, limit, i)] = FALSE
  return(as.numeric(which(p)))
}
d = 3
repeat{
  d = d + 1
  p = getprimes(10^d) # 10^d以下の素数を作る
  p = p[10^(d-1)<p] # d桁の素数だけに絞る
  p_s = do.call(rbind, strsplit(as.character(p), NULL))
  for(n in 0:2){ # n=0,1,2のみ検索
    s = apply(p_s==n, 1, sum) # 各桁の値==nの数
    p_n = p[(p_s[,ncol(p_s)]!=n)&(s>1)&(s%%3==0)] # (一桁目がnではない)かつ(sが3の倍数(0は除く))
    f = unlist(lapply((n+1):9, gsub, pattern=n, x=p_n)) # (n+1) ~ 9に変換
    f = matrix(f%in%p, nrow=length(p_n)) # 素数になるか検索
    f = apply(f, 1, sum) # 素数になる数
    if(any(f>=7))
      return()
  }
}
ans = min(p_n[f>=7])
print(ans)

Problem 52 「置換倍数」

$125874$を2倍すると$251748$となる. これは元の数$125874$と順番は違うが同じ数を含む.
$2x, 3x, 4x, 5x, 6x$が$x$と同じ数を含むような最小の正整数$x$を求めよ.

6倍したときに桁数が変わってはいけないことから、$d$桁のなかで探す場合、$10^{d-1} \sim 10^{d}/6$の範囲内を探せば良い。
例えば$d=3$のときは$100 \sim 167$の範囲。

上限が見えないため、桁数を増やしながら求める数を探す。

d = 2
repeat{
  d = d + 1
  for(n in 10^(d-1):(10^d/6)){
    s = sort(strsplit(as.character(n), NULL)[[1]])
    for(i in 2:7)
      if(any(s!=sort(strsplit(as.character(n*i), NULL)[[1]])))
        break
    if(i==7)
      return()
  }
}
ans = n
print(n)

Problem 53 「組み合わせ選択」

$12345$から3つ選ぶ選び方は10通りである.
$$123, 124, 125, 134, 135, 145, 234, 235, 245, 345.$$組み合わせでは, 以下の記法を用いてこのことを表す:
$C[5,3] = 10$.
一般に,$r ≤ n$について$C[n,r] = n!/(r!(n-r)!)$である. ここで,$n!=n×(n−1)×...×3×2×1,\,0!=1$と階乗を定義する.
$n=23$になるまで, これらの値が100万を超えることはない: $C[23,10] = 1144066$.
$1 ≤ n ≤ 100$について, 100万を超える$C[n, r]$は何通りあるか?

$C[n, r] = C[n-1, r-1] + C[n-1, r]$の関係を使い、動的に求めていく。

limit = 1000000
d = 100
m = diag(d+1)
m[,1] = 1
for(row in 3:(d+1))
  for(col in 2:(row-1))
    m[row, col] = m[row-1, col-1] + m[row-1, col]
ans = sum(m>=(limit))
print(ans)

Problem 54 「ポーカーハンド」

カードゲームのポーカーでは, 手札は5枚のカードからなりランク付けされている. 役を低い方から高い方へ順に並べると以下である.
役無し(ハイカード): 一番値が大きいカード
・ ワン・ペア: 同じ値のカードが2枚
・ ツー・ペア: 2つの異なる値のペア
・ スリーカード: 同じ値のカードが3枚
・ ストレート: 5枚の連続する値のカード
・ フラッシュ: 全てのカードが同じスート (注: スートとはダイヤ・ハート・クラブ/スペードというカードの絵柄のこと)
・ フルハウス: スリーカードとペア
・ フォーカード: 同じ値のカードが4枚
・ ストレートフラッシュ: ストレートかつフラッシュ
・ ロイヤルフラッシュ: 同じスートの10, J, Q, K, A
ここでカードの値は小さい方から2, 3, 4, 5, 6, 7, 8, 9, 10, J, Q, K, Aである. (訳注：データ中で10は'T'と表される)
もし2人のプレイヤーが同じ役の場合には, 役を構成する中で値が最も大きいカードによってランクが決まる: 例えば, 8のペアは5のペアより強い (下の例1を見よ). それでも同じランクの場合には (例えば, 両者ともQのペアの場合), 一番値が大きいカードによってランクが決まる (下の例4を見よ). 一番値が大きいカードが同じ場合には, 次に値が大きいカードが比べれられ, 以下同様にランクを決定する.
例:

試合	プレイヤー1	プレイヤー2	勝者
1	5H 5C 6S 7S KD 5のペア	2C 3S 8S 8D TD 8のペア	プレイヤー2
2	5D 8C 9S JS AC 役無し, A	2C 5C 7D 8S QH 役無し, Q	プレイヤー1
3	2D 9C AS AH AC Aのスリーカード	3D 6D 7D TD QD ダイヤのフラッシュ	プレイヤー2
4	4D 6S 9H QH QC Qのペア, 9	3D 6D 7H QD QS Qのペア, 7	プレイヤー1
5	2H 2D 4C 4D 4S 4-2のフルハウス	3C 3D 3S 9S 9D 3-9のフルハウス	プレイヤー1

poker.txtには1000個のランダムな手札の組が含まれている. 各行は10枚のカードからなる (スペースで区切られている): 最初の5枚がプレイヤー1の手札であり, 残りの5枚がプレイヤー2の手札である. 以下のことを仮定してよい
全ての手札は正しい (使われない文字が出現しない. 同じカードは繰り返されない)
各プレイヤーの手札は特に決まった順に並んでいるわけではない
各勝負で勝敗は必ず決まる
1000回中プレイヤー1が勝つのは何回か? (訳注 : この問題に置いてA 2 3 4 5というストレートは考えなくてもよい)

ハンドの強さを数値に変換し比較する。
役の強さを整数部分とする。
同じ役のときに必要となる構成要素の強さを小数部分とする。

t = read.csv("https://projecteuler.net/project/resources/p054_poker.txt", header=FALSE, sep=" ")
t = as.matrix(t)
t = rbind(t[, 1:5], t[, 6:10])
suit = substr(t, 2, 2)
value = sub("A", "14", sub("K", "13", sub("Q", "12", sub("J", "11", sub("T", "10", substr(t, 1, 1))))))
value = t(apply(value, 1, function(x) sort(as.numeric(x), decreasing=TRUE)))
s = apply(value, 1, function(x) table(x))
rank = rep(NA, nrow(t))
# フラッシュ
f1 = apply(suit, 1, function(x) length(unique(x)))==1
rank[f1] = 5 + apply(value[f1, ], 1, function(x) x%*%10^c(-2, -4, -6, -8, -10))
# ストレート
f2 = ((value[,1] - value[,5])==4) & (apply(value, 1, function(x) length(unique(x)))==5)
rank[f2] = 4 + value[f2, 1]*10^(-2)
# ストレートフラッシュ
rank[f1&f2] = rank[f1&f2]+4
# ワン・ペア
f = unlist(lapply(s, function(x) (length(x)==4)&(max(x)==2)))
rank[f] = 1 + unlist(lapply(s[f], function(x) c(as.numeric(names(which.max(x))), sort(as.numeric(names(x[x==1])), decreasing=TRUE))%*%10^c(-2, -4, -6, -8)))
# ツー・ペア
f = unlist(lapply(s, function(x) (length(x)==3)&(max(x)==2)))
rank[f] = 2 + unlist(lapply(s[f], function(x) c(sort(as.numeric(names(x[x==2])), decreasing=TRUE), as.numeric(names(x[x==1])))%*%10^c(-2, -4, -6)))
# スリーカード
f = unlist(lapply(s, function(x) (length(x)==3)&(max(x)==3)))
rank[f] = 3 + unlist(lapply(s[f], function(x) c(as.numeric(names(which.max(x))), sort(as.numeric(names(x[x==1])), decreasing=TRUE))%*%10^c(-2, -4, -6)))
# フルハウス
f = unlist(lapply(s, function(x) (length(x)==2)&(max(x)==3)))
rank[f] = 6 + unlist(lapply(s[f], function(x) as.numeric(names(sort(table(x), decreasing=TRUE)))%*%10^c(-2, -4)))
# フォーカード:
f = unlist(lapply(s, function(x) (length(x)==2)&(max(x)==4)))
rank[f] = 7 + unlist(lapply(s[f], function(x) as.numeric(names(sort(table(x), decreasing=TRUE)))%*%10^c(-2, -4)))
# ハイカード
f = is.na(rank)
rank[f] = apply(value[f, ], 1, function(x) x%*%10^c(-2, -4, -6, -8, -10))
rank = matrix(rank, ncol=2)
ans = sum(rank[, 1]>rank[, 2])
print(ans)

Problem 55 「Lychrel数」

47とその反転を足し合わせると, $47 + 74 = 121$となり, 回文数になる.
全ての数が素早く回文数になるわけではない. $349$を考えよう,
1. $349 + 943 = 1292$
2. $1292 + 2921 = 4213$
3. $4213 + 3124 = 7337$
$349$は, 3回の操作を経て回文数になる.
まだ証明はされていないが, $196$のようないくつかの数字は回文数にならないと考えられている. 反転したものを足すという操作を経ても回文数にならないものをLychrel数と呼ぶ. 先のような数の理論的な性質により, またこの問題の目的のために, Lychrel数で無いと証明されていない数はLychrel数だと仮定する.
更に, 10000未満の数については，常に以下のどちらか一方が成り立つと仮定してよい.
・50回未満の操作で回文数になる
・まだ誰も回文数まで到達していない
実際, $10677$が50回以上の操作を必要とする最初の数である: $4668731596684224866951378664$(53回の操作で28桁のこの回文数になる).
驚くべきことに, 回文数かつLychrel数であるものが存在する. 最初の数は$4994$である.
10000未満のLychrel数の個数を答えよ.
注: 2007/04/24にLychrel数の理論的な性質を強調するために文面が修正された.

問題通りそのまま計算する。
数値の反転は、文字列に変換せず数値計算で反転させる。
(warningsがでるが答えは一致するため目を瞑る)

limit = 10000-1
ans = 0
for(n in 1:(limit)){
  r = sum(n%/%10^(0:floor(log10(n)))%%10*10^(floor(log10(n)):0))
  for(i in 1:51){
    n = n + r
    r = sum(n%/%10^(0:floor(log10(n)))%%10*10^(floor(log10(n)):0))
    if(n==r)
      break
  }
  if(i==51)
    ans = ans + 1
}
print(ans)

Problem 56 「もっとべき乗の数字和」

Googol($10^{100}$)は非常に大きな数である: 1の後に0が100個続く. $100^{100}$は想像を絶する. 1の後に0が200回続く. その大きさにも関わらず, 両者とも数字和(桁の和)は1である.
$a, b < 100$について自然数$ab$を考える. 数字和の最大値を答えよ.

オーバーフローが発生するため、Problem 16と同様に1桁ずつ分割したベクトルで計算する。

ans = 0
for(a in 1:99){
  v = 1
  for(b in 1:99){
    v = v * a
    r = 0
    for(j in 1:length(v)){
      v[j] = v[j] + r
      r = v[j]%/%10
      v[j] = v[j]%%10
    }
    while(r > 0){
      v = c(v, r%%10)
      r = r%/%10
    }
    if(ans<sum(v))
      ans = sum(v)
  }
}
print(ans)

Problem 57 「平方根の近似分数」

2の平方根は無限に続く連分数で表すことができる.
$$\sqrt{2} = 1 + 1/(2 + 1/(2 + 1/(2 + ... ))) = 1.414213...$$最初の4回の繰り返しを展開すると以下が得られる.
$1 + 1/2 = 3/2 = 1.5$
$1 + 1/(2 + 1/2) = 7/5 = 1.4$
$1 + 1/(2 + 1/(2 + 1/2)) = 17/12 = 1.41666...$
$1 + 1/(2 + 1/(2 + 1/(2 + 1/2))) = 41/29 = 1.41379...$
次の3つの項は$99/70, 239/169, 577/408$である. 第8項は$1393/985$である. これは分子の桁数が分母の桁数を超える最初の例である.
最初の1000項を考えたとき, 分子の桁数が分母の桁数を超える項はいくつあるか?

$n$回繰り返したときの分子を$a_n$、分母を$b_n$とする。
$a_1=3, b_1=2$である。
$a_{n+1} = a_n + 2b_n$
$b_{n+1} = a_n + b_n$
という漸化式が成り立つ。
この漸化式を使い逐次的に求めていく。

オーバーフローが発生するため、Problem 56と同様に1桁ずつ分割したベクトルで計算する。

add = function(a, b){
  a = c(a, rep(0, max(0, length(b)-length(a))))
  b = c(b, rep(0, max(0, length(a)-length(b))))
  v = a + b
  r = 0
  for(j in 1:length(v)){
    v[j] = v[j] + r
    r = v[j]%/%10
    v[j] = v[j]%%10
  }
  while(r > 0){
    v = c(v, r%%10)
    r = r%/%10
  }
  return(v)
}
ans = 0
a = 3
b = 2
for(i in 1:999){
  aa = a
  a = add(aa, 2*b)
  b = add(aa, b)
  if(length(a)>length(b))
    ans = ans + 1
}
print(ans)

Problem 58 「螺旋素数」

1から始めて, 以下のように反時計回りに数字を並べていくと, 辺の長さが7の渦巻きが形成される.
37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18 05 04 03 12 29
40 19 06 01 02 11 28
41 20 07 08 09 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49
面白いことに, 奇平方数が右下の対角線上に出現する. もっと面白いことには, 対角線上の13個の数字のうち, 8個が素数である. ここで割合は$8/13 ≈ 62\%$である.
渦巻きに新しい層を付け加えよう. すると辺の長さが9の渦巻きが出来る. 以下, この操作を繰り返していく. 対角線上の素数の割合が$10\%$未満に落ちる最初の辺の長さを求めよ.

辺の長さが$l$のとき、
右下の数字は$l^2$であり、
四隅のそれぞれの数値はそこから$l-1$ずつ引いていけば良いので、
左下の数字 : $l^2 - (l-1)$
左上の数字 : $l^2 - 2(l-1)$
右上の数字 : $l^2 - 3(l-1)$
となる。
右下の数値が素数になることはない。

isprime = function(n) return(n==2||n==3||all(n%%(2:sqrt(n))!=0))
l = 3
p = 3
while(p/(2*l-1) > 0.1) {
  l = l + 2
  p = p + sum(unlist(lapply(l^2 - (l-1)*(1:3), isprime)))
}
ans = l
print(ans)

Problem 59 「XOR暗号解読」

(訳者注: 文字コードの説明は適当です) 各文字はそれぞれ一意のコードに割り当てられている. よく使われる標準としてASCII (American Standard Code for Information Interchange) がある. ASCIIでは, 大文字A = 65, アスタリスク (*) = 42, 小文字k = 107というふうに割り当てられている.
モダンな暗号化の方法として, テキストファイルの各バイトをASCIIに変換し, 秘密鍵から計算された値とXORを取るという手法がある. XOR関数の良い点は, 暗号化に用いたのと同じ暗号化鍵でXORを取ると平文を復号できる点である. 65 XOR 42 = 107であり, 107 XOR 42 = 65である.
破られない暗号化のためには, 鍵は平文と同じ長さのランダムなバイト列でなければならない. ユーザーは暗号文と暗号化鍵を別々の場所に保存する必要がある. また, もし一方が失われると, 暗号文を復号することは不可能になる.
悲しいかな, この手法はほとんどのユーザーにとって非現実的である. そこで, 鍵の変わりにパスワードを用いる手法が用いられる. パスワードが平文より短ければ (よくあることだが), パスワードは鍵として繰り返し用いられる. この手法では, 安全性を保つために十分長いパスワードを用いる必要があるが, 記憶するためにはある程度短くないといけない.
この問題での課題は簡単になっている. 暗号化鍵は3文字の小文字である. cipher1.txtは暗号化されたASCIIのコードを含んでいる. また, 平文はよく用いられる英単語を含んでいる. この暗号文を復号し, 平文のASCIIでの値の和を求めよ.

問題文の意味が良くわからず、他の方のコードを読みなんとなく分かった。
暗号化鍵は3文字の小文字であるので、26*26*26の全パターンを検索していく。
平文はよく用いられる英単語を含んでいるという意味は、" the "を含んでいると解釈する。
普段使わない関数のを知ることができた。

text_raw = scan("https://projecteuler.net/project/resources/p059_cipher.txt", sep=",")
for(l1 in letters) {
  for(l2 in letters) {
    for(l3 in letters) {
      key = unlist(lapply(c(l1, l2, l3), function(x) strtoi(charToRaw(x), 16L)))
      ascii = bitwXor(text_raw, key)
      text = rawToChar(as.raw(ascii))
      if (grepl(' the ', text))
        return()
    }
  }    
}
ans = sum(ascii)
print(ans)

Problem 60 「素数ペア集合」

素数$3, 7, 109, 673$は非凡な性質を持っている. 任意の2つの素数を任意の順で繋げると, また素数になっている. 例えば, $7$と$109$を用いると, $7109$と$1097$の両方が素数である. これら4つの素数の和は$792$である. これは, このような性質をもつ4つの素数の集合の和の中で最小である.
任意の2つの素数を繋げたときに別の素数が生成される, 5つの素数の集合の和の中で最小のものを求めよ.

最小のものという条件が難しい。難問。
例えば700以下の素数を列挙し、その中から最小意外の条件を満たす4つの素数の集合に$3, 7, 109, 673$があるが、当然$701$や$709$の素数を含む集合で和がさらに小さくなる可能性は0ではない。
しかし、$3, 7, 109, 673$の和である$792$を超える素数を使った集合では、その和は$792$を超えるため最小ではなくなる。

以上のことから、
Step1. 最小の条件を無視してとりあえず条件を満たす5つの素数の集合を探す
Step2. その和以下の素数の中から、条件を満たす5つの素数の集合のうち最小のものを求める
という指針で求める集合を探す。

プログラム中下から2行目のans = solve(ans, set_len, ans)でStep2を行っている。
この行を飛ばすと10数秒で計算でき、同じ答えが求まる。
しかし、その答えが最小になっているという根拠は私は示せず、答えが一致するのは単なる偶然であると思っているため、やはりStep2は行いたい。

getprimes = function(limit){
  p = rep(TRUE, limit)
  p[1] = FALSE
  for(i in c(2, seq(3, sqrt(limit), 2)))
    if(p[i])
      p[seq(i*2, limit, i)] = FALSE
  return(as.numeric(which(p)))
}
solve = function(limit, set_len, ans_limit=limit*set_len){
  p = getprimes(limit)[-c(1, 3)] # 2と5は除外
  p = expand.grid(p1=p, p2=p) # 組み合わせ
  p = p[p$p1<p$p2, ] # 逆順を除外
  p["c1"] = p$p2*10^(floor(log10(p$p1))+1)+p$p1
  p["c2"] = p$p1*10^(floor(log10(p$p2))+1)+p$p2
  p = p[(p$p1+p$p2)<ans_limit, ]
  for(d in getprimes(sqrt(max(p$c1, p$c2)))){ # 結合した値が素数のものだけに絞る
    p = p[(p$c1%%d!=0)|(p$c1==d), ]
    p = p[(p$c2%%d!=0)|(p$c2==d), ]
  }
  p_list = list()
  for(i in unique(p$p1))
    p_list[[as.character(i)]] = as.character(p[p$p1==i, "p2"])
  search = function(p_list, set, candidates, set_len, ans) {
    if(length(set)==set_len)
      return(min(sum(as.numeric(set)), ans))
    for(can in candidates)
      ans = search(p_list, c(set, can), intersect(candidates, p_list[[can]]), set_len, ans)
    return(ans)
  }
  ans = search(p_list, c(), names(p_list), set_len, ans_limit)
  return(ans)
}
limit = 100
set_len = 5
repeat{
  ans = solve(limit, set_len)
  if(ans==limit*set_len)
    limit = limit * 10
  else
    break
}
ans = solve(ans, set_len, ans)
print(ans)