プロジェクトオイラー
http://projecteuler.net/index.php
Q195.
辺の長さが整数で、1つだけ角度が60度で、内接円がr以下の三角形の個数をT(r)とする。T(1053779)を求めよ。
角度60度の頂点の対辺をc、他の辺をa, b(a > b)とすると、余弦定理から、
c2 = a2 - ab + b2
内接円の半径rは、
r = 2S / s = √3ab / 2(a + b + c)
内接円の制限を満たしつつ、上のディオファントス方程式を満たす整数解の数を求める。
上の方程式は、
x2 - xy + y2 = 1
の有理数解を考える。この曲線は必ず(1, 1)を通るから、この点を通って傾きが有理数の直線とこの曲線の交点を考えればよい。直線：y = -k(x - 1) + 1との交点は、
( (k2 + 2k) / (k2 + k + 1), 2k / (k2 + k + 1) )
k = m / nとおいて、分母を払うと、
(a, b, c) = (l(m + 2n)m, l(n + 2m)n, l(m2 + mn + n2))
となる。また、
r = √3lmn / 2
となる。
mとnが素でl=1でも、(a, b, c)が互いに素でないことも注意して数え上げる。
from math import sqrt
import fractionsdef T(r):
s = 0
limit = int(r * sqrt(3) * 2)
for m in xrange(2, limit + 1):
for n in xrange(1, m):
if fractions.gcd(m, n) != 1:
continue
a = (m + 2 * n) * m
b = (n + 2 * m) * n
c = m * m + m * n + n * n
d = fractions.gcd(a, fractions.gcd(b, c))
mn = m * n
if r * sqrt(3) * 2 < mn:
break
s += int(r * d * 2 / (sqrt(3) * mn))
return sprint T(1053779)
