RでProject Euler (Problem 51~60)

[Problem 41~50](https://qiita.com/Rigel_/items/484196dbf1dfa2e391c6)

# 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倍ずつして求める数を探す。

```R
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$の範囲。

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

```R
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]$の関係を使い、動的に求めていく。

```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<br>5のペア|2C 3S 8S 8D TD<br>8のペア|プレイヤー2|
|2|5D 8C 9S JS AC<br>役無し, A|2C 5C 7D 8S QH<br>役無し, Q|プレイヤー1|
|3|2D 9C AS AH AC<br>Aのスリーカード|3D 6D 7D TD QD<br>ダイヤのフラッシュ|プレイヤー2|
|4|4D 6S 9H QH QC<br>Qのペア, 9|3D 6D 7H QD QS<br>Qのペア, 7|プレイヤー1|
|5|2H 2D 4C 4D 4S<br>4-2のフルハウス|3C 3D 3S 9S 9D<br>3-9のフルハウス|プレイヤー1|
[poker.txt](https://projecteuler.net/project/resources/p054_poker.txt)には1000個のランダムな手札の組が含まれている. 各行は10枚のカードからなる (スペースで区切られている): 最初の5枚がプレイヤー1の手札であり, 残りの5枚がプレイヤー2の手札である. 以下のことを仮定してよい
全ての手札は正しい (使われない文字が出現しない. 同じカードは繰り返されない)
各プレイヤーの手札は特に決まった順に並んでいるわけではない
各勝負で勝敗は必ず決まる
1000回中プレイヤー1が勝つのは何回か? (訳注 : この問題に置いてA 2 3 4 5というストレートは考えなくてもよい)

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

```R
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がでるが答えは一致するため目を瞑る)

```R
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桁ずつ分割したベクトルで計算する。

```R
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桁ずつ分割したベクトルで計算する。

```R
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の渦巻きが形成される.
<font color="Red">37</font> 36 35 34 33 32 <font color="Red">31</font>
38 <font color="Red">17</font> 16 15 14 <font color="Red">13</font> 30
39 18 <font color="Red">05</font> 04 <font color="Red">03</font> 12 29
40 19 06 <font color="Red">01</font> 02 11 28
41 20 <font color="Red">07</font> 08 <font color="Red">09</font> 10 27
42 <font color="Red">21</font> 22 23 24 <font color="Red">25</font> 26
<font color="Red">43</font> 44 45 46 47 48 <font color="Red">49</font>
面白いことに, 奇平方数が右下の対角線上に出現する. もっと面白いことには, 対角線上の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)$
となる。
右下の数値が素数になることはない。

```R
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](https://projecteuler.net/project/resources/p059_cipher.txt)は暗号化されたASCIIのコードを含んでいる. また, 平文はよく用いられる英単語を含んでいる. この暗号文を復号し, 平文のASCIIでの値の和を求めよ.

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

```R
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は行いたい。

```R
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)
```