Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
130 lines (100 sloc) 3.05 KB
#
# Solution to Project Euler problem 66
# Copyright (c) Project Nayuki. All rights reserved.
#
# https://www.nayuki.io/page/project-euler-solutions
# https://github.com/nayuki/Project-Euler-solutions
#
import eulerlib, fractions
# Based on this insane theorem: Suppose D > 1 is an integer, non-perfect-square.
#
# Express sqrt(D) as the continued fraction (a0, a1, ..., a_{n-1}, (b0, b1, ..., b_{m-1})),
# where the sequence of b's is the periodic part.
#
# Let p/q (in lowest terms) = (a0, a1, ..., a_{n-1}, b0, b1, ..., b_{m-2}).
# (This is a truncation of the continued fraction with only one period minus the last term.)
#
# Then the minimum solution (x, y) for Pell's equation is given by:
# - (p, q) if m is even
# - (p^2 + D q^2, 2pq) if m is odd
def compute():
ans = max((n for n in range(2, 1001) if (not eulerlib.is_square(n))),
key=smallest_solution_x)
return str(ans)
# Returns the smallest x such that x > 0 and there exists some y such that x^2 - n y^2 = 1.
# Requires n to not be a perfect square.
def smallest_solution_x(n):
contfrac = sqrt_to_continued_fraction(n)
temp = contfrac[0] + contfrac[1][ : -1]
val = fractions.Fraction(temp[-1], 1)
for term in reversed(temp[ : -1]):
val = 1 / val + term
if len(contfrac[1]) % 2 == 0:
return val.numerator
else:
return val.numerator**2 + val.denominator**2 * n
# Returns the periodic continued fraction of sqrt(n). Requires n to not be a perfect square.
# result[0] is the minimal non-periodic prefix, and result[1] is the minimal periodic tail.
def sqrt_to_continued_fraction(n):
terms = []
seen = {}
val = QuadraticSurd(0, 1, 1, n)
while True:
seen[val] = len(seen)
flr = val.floor()
terms.append(flr)
val = (val - QuadraticSurd(flr, 0, 1, val.d)).reciprocal()
if val in seen:
break
split = seen[val]
return (terms[ : split], terms[split : ])
# Represents (a + b * sqrt(d)) / c. d must not be a perfect square.
class QuadraticSurd(object):
def __init__(self, a, b, c, d):
if c == 0:
raise ValueError()
# Simplify
if c < 0:
a = -a
b = -b
c = -c
gcd = fractions.gcd(fractions.gcd(a, b), c)
if gcd != 1:
a //= gcd
b //= gcd
c //= gcd
self.a = a
self.b = b
self.c = c
self.d = d
def __sub__(self, other):
if self.d != other.d:
raise ValueError()
return QuadraticSurd(
self.a * other.c - other.a * self.c,
self.b * other.c - other.b * self.c,
self.c * other.c,
self.d)
def reciprocal(self):
return QuadraticSurd(
-self.a * self.c,
self.b * self.c,
self.b * self.b * self.d - self.a * self.a,
self.d)
def floor(self):
temp = eulerlib.sqrt(self.b * self.b * self.d)
if self.b < 0:
temp = -(temp + 1)
temp += self.a
if temp < 0:
temp -= self.c - 1
return temp // self.c
def __eq__(self, other):
return self.a == other.a and self.b == other.b \
and self.c == other.c and self.d == other.d
def __ne__(self, other):
return not (self == other)
def __hash__(self):
return hash(self.a) + hash(self.b) + hash(self.c) + hash(self.d)
if __name__ == "__main__":
print(compute())