RでProject Euler (Problem 21~30)

[Problem 11~20](https://qiita.com/Rigel_/items/18a5400be1ccdb872f95)
[Problem 31~40](https://qiita.com/Rigel_/items/96ef68f4315a6ec51aa0)

# Problem 21 「友愛数」

>$d(n)$を$n$の真の約数の和と定義する. (真の約数とは$n$以外の約数のことである. )
もし,$d(a) = b$かつ$d(b) = a$($a≠b$のとき) を満たすとき,$a$と$b$は友愛数(親和数)であるという.
例えば,$220$の約数は$1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110$なので$d(220) = 284$である.
また,$284$の約数は$1, 2, 4, 71, 142 なので d(284) = 220$である.
それでは$10000$未満の友愛数の和を求めよ.

エラトステネスの篩と同じ考え方で、$10000$未満の各数に対して真の約数の和を求める。
簡単のため12以下で説明する。
はじめに、各要素が$0$の長さが$12$のベクトルを用意する。
![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/171076/07b9d711-bf45-23cd-375d-7a07b69ffe33.png)
1より大きい1の倍数番目の要素を$+1$する
![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/171076/f76f8736-f495-07da-8e61-a8f758bea736.png)
2より大きい2の倍数番目の要素を$+2$する
![fig3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/171076/54184443-8758-4c89-2ba7-2400e3fa8b65.png)
3より大きい3の倍数番目の要素を$+3$する
![fig4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/171076/89ee1741-04e8-5d7a-98a8-e2b4b3e6f315.png)
同様の計算を$6$までする。
![fig5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/171076/867f41c3-7774-9d5b-e9f9-f5caf1872581.png)
以上の計算で、$12$以下の各数に対して真の約数の和が求まる。

あとは、$d(a) = b$かつ$d(b) = a$を満たすものを足していく。
完全数や同じ友愛数を数えないように注意する。

```R
limit = 10000-1
v = rep(0, limit)
for(i in 1:(limit/2))
  v[seq(2*i, limit, i)] = v[seq(2*i, limit, i)] + i
ans = 0
for(n in 2:limit){
  a = v[n]
  if((a>n)&(a<=limit)){
    b = v[a]
    if(n==b)
      ans = ans + a + b
  }
}
print(ans)
```

# Problem 22 「名前のスコア」

>5000個以上の名前が書かれている46Kのテキストファイル [filenames.txt](https://projecteuler.net/project/resources/p022_names.txt)を用いる. まずアルファベット順にソートせよ.
のち, 各名前についてアルファベットに値を割り振り, リスト中の出現順の数と掛け合わせることで, 名前のスコアを計算する.
たとえば, リストがアルファベット順にソートされているとすると, COLINはリストの938番目にある. またCOLINは 3 + 15 + 12 + 9 + 14 = 53 という値を持つ. よってCOLINは 938 × 53 = 49714 というスコアを持つ.
ファイル中の全名前のスコアの合計を求めよ.

rank関数を使ってランキングを得る。
txtの中には"NA"という名前があり、scan関数の引数na.stringsがデフォルトのままだと、この"NA"さんは欠損として読み込まれる。
そのため、na.stringsに適当な文字列を入れておき、欠損として読み込まれるのを回避する。

```R
t = scan("https://projecteuler.net/project/resources/p022_names.txt", sep=",", quote='"', what="character", na.strings="hogehoge")
r = rank(t)
t = strsplit(t, "")
for(i in 1:length(LETTERS))
  t = lapply(t, gsub, pattern=LETTERS[i], replacement=as.character(i))
t = unlist(lapply(t, function(x){sum(as.numeric(x))}))
ans = sum(t*r)
print(ans)
```

# Problem 23 「非過剰数和」

>完全数とは, その数の真の約数の和がそれ自身と一致する数のことである. たとえば,$28$の真の約数の和は,$1 + 2 + 4 + 7 + 14 = 28$であるので, $28$は完全数である.
真の約数の和がその数よりも少ないものを不足数といい, 真の約数の和がその数よりも大きいものを過剰数と呼ぶ.
$12$は,$1 + 2 + 3 + 4 + 6 = 16$となるので, 最小の過剰数である. よって2つの過剰数の和で書ける最少の数は$24$である. 数学的な解析により, $28123$より大きい任意の整数は2つの過剰数の和で書けることが知られている. 2つの過剰数の和で表せない最大の数がこの上限よりも小さいことは分かっているのだが, この上限を減らすことが出来ていない.
2つの過剰数の和で書き表せない正の整数の総和を求めよ.

Problem 21と同じプログラムで、$28123$以下の真の約数の和を求め、過剰数を探す。
あとは総当りで2つの過剰数の和で表せるか確認していく。

```R
limit = 28123
v = rep(0, limit)
for(i in 1:(limit/2))
  v[seq(2*i, limit, i)] = v[seq(2*i, limit, i)] + i
v = which(v>(1:limit))
t = rep(FALSE, limit*2)
for(i in 1:length(v))
  t[v[i]+v[1:i]] = TRUE
ans = sum(which(!head(t, limit)))
print(ans)
```

# Problem 24 「辞書式順列」

>順列とはモノの順番付きの並びのことである. たとえば,$3124$は数$1, 2, 3, 4$の一つの順列である. すべての順列を数の大小でまたは辞書式に並べたものを辞書順と呼ぶ. $0$と$1$と$2$の順列を辞書順に並べると
$012$ $021$ $102$ $120$ $201$ $210$
になる.
$0,1,2,3,4,5,6,7,8,9$からなる順列を辞書式に並べたときの100万番目はいくつか?

例えば先頭に$0$がつくものは$9!=362880$個ある。
従って、先頭に$2$がつくものは$362880×2+1$番目から$362880×3$番目なので100万番目はこの中に入るので、先頭は$2$である。
この考え方で、上位の桁から一桁ずつ決定していく。

```R
n = as.character(0:9)
m = 1000000
ans = c()
for(i in 1:length(n)){
  p = gamma(length(n))
  a = ceiling(m/p)
  ans = c(ans, n[a])
  n = n[-a]
  m = m - p * (a-1)
}
ans = paste(ans, collapse="")
print(ans)
```

# Problem 25 「1000桁のフィボナッチ数」

>フィボナッチ数列は以下の漸化式で定義される:
$F_n = F_{n-1} + F_{n-2}$, ただし$F_1 = 1, F_2 = 1$.
最初の12項は以下である.
$F_1 = 1$
$F_2 = 1$
$F_3 = 2$
$F_4 = 3$
$F_5 = 5$
$F_6 = 8$
$F_7 = 13$
$F_8 = 21$
$F_9 = 34$
$F_{10} = 55$
$F_{11} = 89$
$F_{12} = 144$
12番目の項, $F_{12}$が3桁になる最初の項である.
1000桁になる最初の項の番号を答えよ.

Problem 2の考え方で一発で決定する。

```R
l = 1000
phi = (1+sqrt(5))/2
ans = ceiling((l-1 + log10(sqrt(5)))/log10(phi))
print(ans)
```

# Problem 26 「逆数の循環節 その1」

>単位分数とは分子が1の分数である. 分母が2から10の単位分数を10進数で表記すると次のようになる.
$1/2 = 0.5$
$1/3 = 0.(3)$
$1/4 = 0.25$
$1/5 = 0.2$
$1/6 = 0.1(6)$
$1/7 = 0.(142857)$
$1/8 = 0.125$
$1/9 = 0.(1)$
$1/10 = 0.1$
$0.1(6)$は$0.166666...$という数字であり, 1桁の循環節を持つ.$1/7$の循環節は6桁ある.
$d < 1000$なる$1/d$の中で小数部の循環節が最も長くなるような$d$を求めよ.

$1/7$について考える。以下のような繰り返し計算を行う。
　1回目：$10\equiv 3\ \pmod{7}$
　2回目：$3*10\equiv 2\ \pmod{7}$
　3回目：$2*10\equiv 6\ \pmod{7}$
　4回目：$6*10\equiv 4\ \pmod{7}$
　5回目：$4*10\equiv 5\ \pmod{7}$
　6回目：$5*10\equiv 1\ \pmod{7}$
ここで、あまりが1となったので、7回目は1回目と同じ計算になり、以後繰り返される。
この循環節の長さは6桁になる。
このような計算を1000以下の素数で行っていく。
参考→[循環小数と循環節の長さ](http://www.virtual-hs.com/math/seisu009.html)

```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)))
}
limit = 1000
p = getprimes(limit)
r = rep(1, length(p))
ans = c()
for(i in 1:limit){
  r = r * 10
  r = r %% p
  f = r == 1
  pp = p[f]
  p = p[!f]
  r = r[!f]
  ans = c(ans, pp[pp==i+1])
}
ans = max(ans)
print(ans)
```

# Problem 27 「二次式素数」

>オイラーは以下の二次式を考案している:$$n^2 + n + 41.$$この式は,$n$を$0$から$39$までの連続する整数としたときに40個の素数を生成する. しかし,$n = 40$のとき$40^2 + 40 + 41=40(40 + 1) + 41$となり$41$で割り切れる. また,$n = 41$のときは$41^2 + 41 + 41$であり明らかに$41$で割り切れる.
計算機を用いて, 二次式$n^2 - 79n + 1601$という式が発見できた. これは$n = 0$から$79$の連続する整数で80個の素数を生成する. 係数の積は,$-79 × 1601$で$-126479$である.
さて,$|a| < 1000, |b| ≤ 1000$として以下の二次式を考える (ここで$|a|$は絶対値): 例えば$|11| = 11 |-4| = 4$である.
$$n^2 + an + b$$$n = 0$から始めて連続する整数で素数を生成したときに最長の長さとなる上の二次式の, 係数$a, b$の積を答えよ.

$n=0$のときを考えると、$b$は素数。
$n=2, b=2$のとき、$n^2 + an + b$は偶数となるので$b$は$2$ではない。
２つの条件を合わせると、$b$は$3$以上の素数。
$a$が偶数、$n$が奇数のとき、$n^2 + an + b$は偶数となるので$a$は奇数。
$n=1$のとき、$n^2 + an + b = 1 + a + b \geq 3$より、$2-b\leq a \leq999$。
以上の条件で絞れたパターンをすべて調べる。


```R
limit = 2*1000^2
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
m = 0
for(b in which(head(p, 1000))[-1]){
  for(a in seq(2-b, 999, by=2)){
    for(n in 1:(b-1)){
      if(n^2+a*n+b<2) break
      if(!p[n^2+a*n+b]) break
    }
    if(m<=n){
      m = n
      ans = a*b
    }
  }
}
print(ans)
```

# Problem 28 「螺旋状に並んだ数の対角線」

>1から初めて右方向に進み時計回りに数字を増やしていき, 5×5の螺旋が以下のように生成される:
<font color="Red">21</font> 22 23 24 <font color="Red">25</font>
20 <font color="Red">07</font> 08 <font color="Red">09</font> 10
19 06 <font color="Red">01</font> 02 11
18 <font color="Red">05</font> 04 <font color="Red">03</font> 12
<font color="Red">17</font> 16 15 14 <font color="Red">13</font>
両対角線上の数字の合計は101であることが確かめられる.
1001×1001の螺旋を同じ方法で生成したとき, 対角線上の数字の和はいくつか?

$n\geq 3$のとき、$n×n$の、
右上の数字は：$n^2$
左上の数字は：$n^2-(n-1)$
左下の数字は：$n^2-2(n-1)$
右下の数字は：$n^2-3(n-1)$
よって、その4つの数字の和は$$n^2 + n^2-(n-1) + n^2-2(n-1) + n^2-3(n-1) = 4n^2 - 6n + 6$$である。
これを$n=3, 5, .., 1001$で足し上げる。

```R
ans = 1
for(n in seq(3, 1001, by=2))
  ans = ans + 4*n^2 - 6*n + 6
print(ans)
```

# Problem 29 「個別のべき乗」

>$2 ≤ a ≤ 5$と$2 ≤ b ≤ 5$について,$a^b$を全て考えてみよう:
$2^2=4, 2^3=8, 2^4=16, 2^5=32$
$3^2=9, 3^3=27, 3^4=81, 3^5=243$
$4^2=16, 4^3=64, 4^4=256, 4^5=1024$
$5^2=25, 5^3=125, 5^4=625, 5^5=3125$
これらを小さい順に並べ, 同じ数を除いたとすると, 15個の項を得る:
$4, 8, 9, 16, 25, 27, 32, 64, 81, 125, 243, 256, 625, 1024, 3125$
$2 ≤ a ≤ 100, 2 ≤ b ≤ 100$で同じことをしたときいくつの異なる項が存在するか?

$100×100=10000$パターンくらいしかないので、全パターン生成し、uniqueとって長さを求める。

```R
d = expand.grid(2:100, 2:100)
ans = length(unique(d$Var1^d$Var2))
print(ans)
```

# Problem 30 「各桁の5乗」

>驚くべきことに, 各桁を4乗した数の和が元の数と一致する数は3つしかない.
$1634 = 1^4 + 6^4 + 3^4 + 4^4$
$8208 = 8^4 + 2^4 + 0^4 + 8^4$
$9474 = 9^4 + 4^4 + 7^4 + 4^4$
ただし, $1=1^4$は含まないものとする. この数たちの和は$1634 + 8208 + 9474 = 19316$である.
各桁を5乗した数の和が元の数と一致するような数の総和を求めよ.

$7*9^5 = 413343$：6桁
$6*9^5 = 354294$：6桁
なので、求めたい数は6桁以下である。
特に工夫もなく$0$から$999999$まで全パターンから条件に一致する数を探す。

```R
v = 0:(10^6-1)
d = Reduce("+", lapply(1:6, function(x) rep(rep((0:9)^5, each=10^(x-1)), times=10^(6-x))))
ans = sum(v[v == d])-1
print(ans)
```