プロジェクトオイラー
http://projecteuler.net/index.php
Q66.
ペル方程式 x2 - Dy2 = 1で最小の解xが最も大きいD≤1000
ここを参照した。この計算をするために、前に作った代数的数のクラスを少し拡張した。
# coding: s_jisfrom math import sqrt
class Gal:
def __init__(self, a, b = 0, c = 0, d = 1):
self.a, self.b, self.c, self.d = a, b, c, d
self.normalize()
def __add__(self, a):
if isinstance(a, Gal):
f = self.GCD(self.d, a.d)
f1 = a.d / f
f2 = self.d / f
return Gal(f1 * self.a + f2 * a.a,
f1 * self.b + f2 * a.b, self.c, f2 * a.d)
else: # 第2引数は整数とみなす
return Gal(self.a + a * self.d, self.b, self.c, self.d)
def __radd__(self, a):
return Gal(self.a + a * self.d, self.b, self.c, self.d)
def __sub__(self, a):
if isinstance(a, Gal):
f = self.GCD(self.d, a.d)
f1 = a.d / f
f2 = self.d / f
return Gal(f1 * self.a - f2 * a.a,
f1 * self.b - f2 * a.b, self.c, f2 * a.d)
else: # 第2引数は整数とみなす
return Gal(self.a - a * self.d, self.b, self.c, self.d)
def __rsub__(self, a):
return Gal(-self.a + a * self.d, -self.b, self.c, self.d)
def __mul__(self, a):
if isinstance(a, Gal):
return Gal(self.a * a.a + self.b * a.b * self.c,
self.a * a.b + self.b * a.a, self.c, self.d * a.d)
else:
return Gal(self.a * a, self.b * a, self.c, self.d)
def __rmul__(self, a):
return Gal(self.a * a, self.b * a, self.c, self.d)
def __div__(self, a):
if isinstance(a, Gal):
square = a.a * a.a - a.b * a.b * a.c
return Gal((self.a * a.a - self.b * a.b * self.c) * a.d,
(-self.a * a.b + self.b * a.a) * a.d,
self.c, self.d * square)
else:
return Gal(self.a, self.b, self.c, self.d * a)
def __rdiv__(self, a):
square = self.a * self.a - self.b * self.b * self.c
a *= self.d
return Gal(a * self.a, -a * self.b, self.c, square)
def __neg__(self):
return Gal(-self.a, -self.b, self.c, self.d)
def conjugate(self):
return Gal(self.a, -self.b, self.c, self.d)
def val(self):
return (self.a + self.b * sqrt(self.c)) / self.d
def normalize(self):
f = self.GCD(self.GCD(self.a, self.b), self.d)
if self.d < 0:
f = -f
self.a /= f
self.b /= f
self.d /= f
def __str__(self):
if self.d == 1:
return self.str_core()
elif self.a == 0 or self.b == 0:
return self.str_core() + "/" + str(self.d)
else:
return "(" + self.str_core() + ")/" + str(self.d)
def str_core(self):
if self.a == 0 and self.b == 0:
result = "0",
else:
result = ""
if self.a != 0:
result = str(self.a)
if self.a != 0 and self.b > 0:
result += "+"
if self.b == 1:
result += "√" + str(self.c)
elif self.b == -1:
result += "-√" + str(self.c)
elif self.b != 0:
result += str(self.b) + "√" + str(self.c)
return result;
def GCD(self, a, b):
if a == 0:
return b
elif b == 0:
return a
a1 = abs(a)
b1 = abs(b)
while True:
a2 = a1 % b1
if a2 == 0:
return b1
b1 %= a1
if b1 == 0:
return a1
a1 = a2
from itertools import count, takewhile
from math import sqrt
from gal import Galdef is_square(n):
m = int(sqrt(n + 0.5))
return m * m == ndef gen_not_square():
for n in count(2):
if not is_square(n):
yield ndef min_solution(D):
a0 = Gal(0, 1, D, 1)
a1 = Gal(1, 0, D)
for n in count():
b = -a0.conjugate() / a1.conjugate()
if n > 1 and a1.a ** 2 - D * a1.b ** 2 == 1:
return a1.a
first = False
k = int(b.val())
a0, a1 = a1, a0 + k * a1print reduce(lambda x, y: max(x, y, key = lambda x: x[1]),
map(lambda D: [ D, min_solution(D) ],
takewhile(lambda n: n <= 1000, gen_not_square())))[0]