Problem 136 - もうカツ丼はいいよな

Problem 136 - Project Euler

正整数かつ等差数列である3つの数x, y, zと、正整数nについての方程式 $x^2-y^2-z^2=n$ を考える。n=20のとき、方程式はただ一つの解を持つ。
$13^2 - 10^2 - 7^2 = 20$
100未満のnについて、そのうち25個が方程式に対しだた一つの解を持つ。
5千万より小さいnについて、ただ一つの解を持つものがいくつあるか答えよ。

解答

nの値としてあり得るのは、pを奇素数または1として

$n = p$ ( $n \equiv 3 (\bmod 4)$ )
$n = 4p$
$n = 16p$

の3パターンしかない。なので、5千万までの素数表さえ用意出来ればすぐに解ける。フォーラム覗いてみてもこのやり方の人が多い。といっても5千万以下ともなるとパッと求めるのは厳しいのでなんか別の方法があるかも。

証明

等差数列の公差をaとおき、整理すると
$\frac{n}{y}+y=4a$
また $y \geq a+1$ , すなわち $\frac{n}{y}+y \leq 4y-4$ 。
よって
$\frac{n}{y}+y \equiv 0 (\bmod 4)$
となり、上記不等式も満たすようなyがただ一つ存在するようなnの条件を求めればよい。

nが素数のとき

$n/y$ は整数でなければならないので、yとしては1またはnのみがあり得る。y=1では不等式を満たさないが、y=nならば必ず満たす。よって
$1+n=4a$
すなわち
$n = p$ ( $n \equiv 3 (\bmod 4)$ )
を満たすnは全てただ一つの解を持つ。

nが合成数のとき

nが異なる2つの奇数o, rに分解できたとして、y=oのとき解が出来たとする。
$\frac{or}{o} + o = 4a$
整理すると
$r = 4a - o$
ここで、y = n = orとした上、上式を代入する。
$1 + 4oa - o^2 = 4b$
pは奇数なので、 $2m-1$ の様に書ける。これを代入すると
$4pa - 4(m^2 - m)$
が得られるので、y=nも解を作る。よって解は2つ以上となる。
以上より、nは約数として2つ以上の奇数を含まず、次の形に素因数分解できる。
$n = 2^ip$
よって
$\frac{n}{y}+y=4a$
は
$\frac{2^ip}{2^j} + 2^j = 4a$
もしくは
$\frac{2^ip}{2^jp} + 2^jp = 4a$
と書ける。この式を成立させるためには、 $i=2$ かつ $j=1$ もしくは $i-j \geq 2$ かつ $j \geq 2$ のいずれかでなければならない。
不等式を考慮すると、 $i=2$ かつ $j=1$ のときは $\frac{2^ip}{2^jp} + 2^jp = 4a$ のみが必ず成立する。つまり、 $n=4p$ は必ず解を一つ作る。
また、 $i \geq 5$ では2つ以上の解を必ず作ることができる。
よって $i=4$ のみが残る可能性となる。
yとしてあり得るのは10通りしかないのですべて調べると、 $y=4p$ のときに必ず、またそのときのみ解が出来ることがわかる。よって $n=16p$ も必ず解を一つ作る。(証明終)

ソース

limit <- 5e7
sieve <- logical(limit)
sieve[1] <- TRUE
sieve[seq(4, limit, by=2)] <- TRUE
for(i in seq(3, sqrt(limit), by=2)){
  sieve[seq(i*2, limit, i)] <- TRUE
}
prime <- c(1, (1:limit)[!sieve][-1])

sum(prime %% 4 == 3) + 
  sum(4 * prime < limit) + 
  sum(16 * prime < limit)

45秒程度。ほとんどは素数表を求める時間で、残りは1sec以下で計算できる。
(追記)エラトステネスの篩を見直したら6.7秒に短縮できた。

limit <- 5e7

sieve <- logical((limit-1)/2)
for(i in 1:(floor(sqrt(limit)-1)/2)){
  if(!sieve[i]){
    sieve[seq(2*i*(i+1), length(sieve), by = 2*i+1)] <- TRUE
  }
}
prime <- c(2, (1:length(sieve))[!sieve] * 2 + 1)

sum(prime %% 4 == 3) + 
  sum(4 * prime < limit) + 
  sum(16 * prime < limit)