プロジェクトオイラー
http://projecteuler.net/index.php
Q198.
9/40のように分母5以下で近似すると、1/4と1/5の2つ出てくることがある。1/100未満で、分母が108を超えないとき、そのような正数はいくつかるか。
恐らく6日間考えた。昨日になってようやく糸口がつかめた。
例のように、1/4と1/5の間には明らかにこれらより分母が小さい分数はない。あるとすれば分子が1でなければならないが、4と5の間に整数はない。よって、これらの平均である9/40の近似が2つ出てくることになる。
もっと一般的に、f1 = p1/q1とf2 = p2/q2の間に分母が小さな分数がない必要十分条件は、
p2q1 - p1q2 = 1
である。
例えば、f1の方を固定すると、一次不定方程式を解かなければならない。そうするとあらかじめ最小解の大雑把な大きさの予想がつかないから、かなり大きな範囲を探索しなければならない。これはQ.251と同じである。そこで、f2が特定の値になるf1を選んで解くことにする。こうすると、結局は一次不定方程式を解かなければならないことは同じなのだが、探索する範囲がグッと減る。
最後に、重複する解を取り除くのに苦しんだが、結局安易にsetを使った。もう少し考えればよかったかもしれないが、30秒で終わったのでよいことにした。
この問題は正答者数からみても、最強の問題の一つだと思う。240番台に正答者が少ない問題があるが、多くの人は恐らくそこまではたどり着いていない。Q.251は未だに63人で、これも最強問題の一つだと思う。
import fractionsdef linear_diophantine(a, b, c):
if a == 1:
return (b + c, 1)
elif a == -1:
return (-b - c, 1)
elif a == 0:
return (1, -c / b)
d = b / a
r = b % a
t = linear_diophantine(r, a, -c)
return (t[1] + d * t[0], t[0])def count(m, n):
s = set()
for q0 in xrange(m, int( (n / 2 + 0.5) ** 0.5)):
for p0 in xrange(q0 / m + 1):
if fractions.gcd(p0, q0) == 1:
q10, p10 = linear_diophantine(p0, q0, 1)
kmax = (n / (2 * q0) - q10) / q0
for k in xrange(kmax + 1):
p1 = p10 + k * p0
q1 = q10 + k * q0
# s += (n - 2 * q0 * q1) / (2 * q1 * q1) + 1
for l in xrange((n - 2 * q0 * q1) / (2 * q1 * q1) + 1):
p2 = p0 + l * p1
q2 = q0 + l * q1
f = (p1 * q2 + p2 * q1, 2 * q1 * q2)
s.add(f)
return len(s)M = 100
N = 10 ** 8
s = N / 2 - M / 2 + count(M, N)
print s
