RでProject Euler (Problem 41~50)

Problem 31~40
Problem 51~60

Problem 41 「パンデジタル素数」

$n$桁パンデジタルであるとは, $1$から$n$までの数を各桁に1つずつ持つこととする.
下のリンク先にあるような数学的定義とは異なる.
例えば$2143$は$4$桁パンデジタル数であり, かつ素数である. $n$桁（この問題の定義では$9$桁以下）パンデジタルな素数の中で最大の数を答えよ.

$8$桁パンデジタル数は、各桁の数の和が$1+2+...+8=36$なので3の倍数となり素数にはならない。
$9$桁パンデジタル数は、各桁の数の和が$1+2+...+9=45$なので3の倍数となり素数にはならない。
なので、$7$桁以下で考える。
$7$桁の最大のパンデジタル数は$7654321$であり、この値以下の素数をエラトステネスの篩で求めておく。
そして、その素数の中からパンデジタル数を選ぶ。

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 = 7654321
p = getprimes(limit)
p = p[!grepl(0, p)]
for(i in 1:7)
  p = p[((log10(p)+1)<i) | grepl(i, p)]
ans = max(p)
print(ans)

Problem 42 「符号化三角数」

三角数の$n$項は$t_n = n(n+1)/2$で与えられる. 最初の$10$項は
$1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...$
である.
単語中のアルファベットを数値に変換した後に和をとる. この和を「単語の値」と呼ぶことにする. 例えば SKY は$19 + 11 + 25 = 55 = t_{10}$である. 単語の値が三角数であるとき, その単語を三角語と呼ぶ.
16Kのテキストファイル words.txt中に約2000語の英単語が記されている. 三角語はいくつあるか?

ある自然数$N$が三角数であるかの判定を考える。
$N = t_n = n(n+1)/2$ を$n$について解くと、
$$n = \frac{\pm \sqrt{8N+1}-1}{2}$$となる。
ここで、プラスマイナスのうち、マイナスはn<0となるため不適。
$n$が整数となるためには、$\sqrt{8N+1}$が整数でなければいけない。
そのとき、$\sqrt{8N+1}-1$は必ず偶数となるため、$n$は必ず整数となる。
従って、$\sqrt{8N+1}$が整数となるかを判定することで、自然数$N$が三角数であるか判定することができる。

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

Problem 43 「部分文字列被整除性」

数$1406357289$は$0$から$9$のパンデジタル数である ($0$から$9$が1度ずつ現れるので). この数は部分文字列が面白い性質を持っている.
$d_1$を上位$1$桁目, $d_2$を上位$2$桁目の数とし, 以下順に$d_n$を定義する. この記法を用いると次のことが分かる.
$d_2d_3d_4=406$は$2$で割り切れる
$d_3d_4d_5=063$は$3$で割り切れる
$d_4d_5d_6=635$は$5$で割り切れる
$d_5d_6d_7=357$は$7$で割り切れる
$d_6d_7d_8=572$は$11$で割り切れる
$d_7d_8d_9=728$は$13$で割り切れる
$d_8d_9d_{10}=289$は$17$で割り切れる
このような性質をもつ$0$から$9$のパンデジタル数の総和を求めよ.

・$3$桁の数を生成
・そのうち下$3$桁が$2$で割り切れるものだけに絞る
・右に$1$桁追加
・そのうち下$3$桁が$3$で割り切れるものだけに絞る
・...
という処理を繰り返し、最後にパンデジタル数のものだけ抽出する。

p = c(2, 3, 5, 7, 11, 13, 17)
v = 123:987
for(i in 1:length(p)){
  v = c(outer(10*v, 0:9, "+"))
  v = v[v%%1000%%p[i]==0]
}
v = v[Reduce("&", lapply(0:9, function(i) grepl(i, v)))]
ans = sum(v)
print(ans)

Problem 44 「五角数」

五角数は$P_n = n(3n-1)/2$で生成される. 最初の$10$項は
$1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ...$
である.
$P_4 + P_7 = 22 + 70 = 92 = P_8$である. しかし差$70 - 22 = 48$は五角数ではない.
五角数のペア$P_j$と$P_k$について, 差と和が五角数になるものを考える. 差を$D = |P_k - P_j|$と書く. 差$D$の最小値を求めよ.

ある自然数$N$が五角数であるかは、$\sqrt{24N+1}+1$が$6$で割り切れるかどうかで判定することができる。
while文で五角数を作っていき、逐次的に判定する。

penta = c()
i = 1
p = 1
f = FALSE
while(!any(f)){
  penta = c(penta, p)
  i = i + 1
  p = i*(3*i-1)/2
  f = ((sqrt(24*(p+penta)+1)+1)%%6==0)&((sqrt(24*(p-penta)+1)+1)%%6==0)
}
ans = p - penta[f]
print(ans)

Problem 45 「三角数, 五角数, 六角数」

三角数, 五角数, 六角数は以下のように生成される.
三角数: $T_n=n(n+1)/2$,　　　$\{1, 3, 6, 10, 15, ...\}$
五角数: $P_n=n(3n-1)/2$,　　　$\{1, 5, 12, 22, 35, ...\}$
六角数: $H_n=n(2n-1)$,　　　$\{1, 6, 15, 28, 45, ...\}$
$T_{285} = P_{165} = H_{143} = 40755$であることが分かる.
次の三角数かつ五角数かつ六角数な数を求めよ.

Problem42とProblem44で用いた三角数と五角数の判定方法を使用する。
while文で六角数を作っていき、逐次的に判定する。

i=144
ans = i*(2*i-1)
while(((sqrt(24*ans+1)+1)%%6!=0)||(sqrt(8*ans+1)%%1!=0)){
  i = i + 1
  ans = i*(2*i-1)
}
print(ans)

Problem 46 「もうひとつのゴールドバッハの予想」

Christian Goldbachは全ての奇合成数は平方数の2倍と素数の和で表せると予想した.
$9 = 7 + 2×1^2$
$15 = 7 + 2×2^2$
$21 = 3 + 2×3^2$
$25 = 7 + 2×3^2$
$27 = 19 + 2×2^2$
$33 = 31 + 2×1^2$
後に, この予想は誤りであることが分かった.
平方数の2倍と素数の和で表せない最小の奇合成数はいくつか?

上限が見えないため、上限を10倍ずつしていき求める数を探す。
上限以下の素数をエラトステネスの篩で作成しておく。
そして、"奇合成数 - 2×平方数"を全パターン作成し、それが素数になるかどうかを調べる。

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)))
}
f = TRUE
limit = 1
while(all(f)){
  limit = limit * 10
  p = getprimes(limit)
  v = seq(3, limit, 2)
  v = v[!v%in%p]
  f = unlist(lapply(v, function(x) any((x-2*(1:sqrt(x/2))^2)%in%p)))
}
ans = min(v[!f])
print(ans)

Problem 47 「異なる素因数」

それぞれ2つの異なる素因数を持つ連続する2つの数が最初に現れるのは:
$14 = 2 × 7$
$15 = 3 × 5$
それぞれ3つの異なる素因数を持つ連続する3つの数が最初に現れるのは:
$644 = 2^2 × 7 × 23$
$645 = 3 × 5 × 43$
$646 = 2 × 17 × 19$
最初に現れるそれぞれ4つの異なる素因数を持つ連続する4つの数を求めよ. その最初の数はいくつか?

ある自然数$n$が異なる素因数をいくつもつかを返す関数を作成する。
あとは愚直にwhile文で4つ連続で4つの異なる素因数を持つ数を探す。

l = 4
count = function(n){
  i = 0
  for(d in c(2, seq(3, sqrt(n), 2))){
    if(n%%d==0)
      i = i + 1
    while(n%%d==0)
      n = n / d
    if(n==1)
      return(i)
  }
  return(i+1)  
}
cnt = 0
n = 9
while(cnt!=l){
  if(count(n)==l)
    cnt = cnt + 1
  else
    cnt = 0
  n = n + 1
}
ans = n-l
print(ans)

Problem 48 「自身のべき乗(self powers)」

次の式は,$1^1 + 2^2 + 3^3 + ... + 10^{10} = 10405071317$である.
では,$1^1 + 2^2 + 3^3 + ... + 1000^{1000}$の最後の10桁を求めよ.

普通に計算するとオーバーフローが発生し、正しい答えを求めることができない。
modの性質を利用して、累乗する過程で都度$10^{10}$で割った余りにしておく。

p = function(n){
  ans = 1
  for(i in 1:n)
    ans = (ans * n)%%(10^10)
  return(ans)
}
ans = sum(unlist(lapply(1:1000, p)))%%(10^10)
print(ans)

Problem 49 「素数数列」

項差$3330$の等差数列$1487, 4817, 8147$は次の2つの変わった性質を持つ.
(i)3つの項はそれぞれ素数である.
(ii)各項は他の項の置換で表される.
1, 2, 3桁の素数にはこのような性質を持った数列は存在しないが, 4桁の増加列にはもう1つ存在する.
それではこの数列の3つの項を連結した12桁の数を求めよ.

はじめに、$1000$以上$10000$以下の素数をエラトステネスの篩で求めておく。
ある$i$番目の素数$p_i$に対して$p_i+d$と$p_i+2d$も素数となる交差$d$は、
・偶数
・$p_{i+1} - p_i \leq d \leq (10000-p_i)/2$
の条件を満たす。
そういった$d$を列挙し等差数列を作成する。
その等差数列のうち"各項は他の項の置換で表される."ものを抽出する。

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 = 10^4
p = getprimes(limit)
p = p[p>1000]
ans = c()
for(i in 1:(length(p)-1)){
  d = seq(p[2]-p[1], (limit-p[1])%/%2, 2)
  d = d[((p[1]+d)%in%p) & ((p[1]+2*d)%in%p)]
  ps = sort(strsplit(as.character(p[1]), "")[[1]])
  f = unlist(lapply(strsplit(as.character(p[1]+d), ""), function(x) all(sort(x)==ps)))
  d = d[f]
  f = unlist(lapply(strsplit(as.character(p[1]+2*d), ""), function(x) all(sort(x)==ps)))
  if(any(f) & p[1]!=1487)
    ans = c(ans, paste(p[1]+c(0, d[f], 2*d[f]), collapse=""))
  p = p[-1]
}
print(ans)

Problem 50 「連続する素数の和」

素数$41$は6つの連続する素数の和として表せる:
$41 = 2 + 3 + 5 + 7 + 11 + 13$.
100未満の素数を連続する素数の和で表したときにこれが最長になる.
同様に, 連続する素数の和で1000未満の素数を表したときに最長になるのは$953$で21項を持つ.
100万未満の素数を連続する素数の和で表したときに最長になるのはどの素数か?

はじめに、100万未満の素数をエラトステネスの篩で求めておく。
あとは効率良く全探索すると意外に早く終る。
現在$i$番目の素数$p_i$までチェックしており、見つけている最長の項数を$L$項とすると、
$p_{i+1}$から$p_{i+L}$の和が100万を超えるとストップして良い。
ちなみに、求める連続する素数は$7$から始まる。
そのため、$2, 3, 5, 7$の確認でループを抜ける。

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 = 10^6
p = getprimes(limit)
L = 1
while(sum(p[1:L])<limit){
  l = max(which(cumsum(p)%in%p))
  if(L<l){
    L=l
    ans = sum(p[1:l])
  }
  p = p[-1]
}
print(ans)