98 lines
2.4 KiB
Python
98 lines
2.4 KiB
Python
import math
|
|
from e057 import add_fractions
|
|
|
|
|
|
def get_floor_sqrt(n):
|
|
return math.floor(math.sqrt(n))
|
|
|
|
|
|
def next_expansion_1(current_a, current_nominator,
|
|
current_denominator, original_number):
|
|
cn = current_nominator
|
|
cd = current_denominator
|
|
cn, cd = cd, (original_number - cd * cd) // cn
|
|
next_a = math.floor((math.sqrt(original_number) + cn) // cd)
|
|
cn = cn - next_a * cd
|
|
cn *= -1
|
|
|
|
return next_a, cn, cd
|
|
|
|
|
|
def get_continued_fraction_sequence(n):
|
|
|
|
# If number is a square number we return it.
|
|
floor_sqrt = get_floor_sqrt(n)
|
|
if n == floor_sqrt * floor_sqrt:
|
|
return ((floor_sqrt), [])
|
|
|
|
# Otherwise, we calculate the next expansion till we
|
|
# encounter a step a second time.
|
|
a = floor_sqrt
|
|
cn = a
|
|
cd = 1
|
|
sequence = []
|
|
previous_steps = []
|
|
|
|
while not (a, cn, cd) in previous_steps:
|
|
# print("a: {} cn: {} cd: {}".format(a, cn, cd))
|
|
previous_steps.append((a, cn, cd))
|
|
a, cn, cd = next_expansion_1(a, cd, cn, n)
|
|
sequence.append(a)
|
|
sequence.pop()
|
|
return ((floor_sqrt), sequence)
|
|
|
|
|
|
def next_expansion_2(previous_numerator, previous_denumerator, value):
|
|
if previous_numerator == 0:
|
|
return (value, 1)
|
|
return add_fractions(previous_denumerator, previous_numerator, value, 1)
|
|
|
|
|
|
def get_fractions(n, x):
|
|
# print("get_fractions(n={}, x={})".format(n, x))
|
|
# Get sequence of a_x
|
|
first_value, sequence = get_continued_fraction_sequence(n)
|
|
sequence = [first_value] + math.ceil(x / len(sequence)) * sequence
|
|
sequence = sequence[:x + 1]
|
|
sequence = sequence[::-1]
|
|
|
|
n, d = 0, 1
|
|
for s in sequence:
|
|
n, d = next_expansion_2(n, d, s)
|
|
return (n, d)
|
|
|
|
|
|
def get_minimal_solution(d):
|
|
for i in range(0, 100):
|
|
x, y = get_fractions(d, i)
|
|
if x * x - d * y * y == 1:
|
|
return((x, y))
|
|
|
|
|
|
def is_square(n):
|
|
return math.sqrt(n).is_integer()
|
|
|
|
|
|
def euler_066():
|
|
# XXX: This seems really long and complicated. We can do better.
|
|
x_max = 0
|
|
d_max = 0
|
|
|
|
for d in range(2, 1001):
|
|
if is_square(d):
|
|
continue
|
|
|
|
x, y = get_minimal_solution(d)
|
|
if x > x_max:
|
|
x_max = x
|
|
d_max = d
|
|
|
|
print("d: {} x: {}".format(d_max, x_max))
|
|
s = d_max
|
|
return s
|
|
|
|
|
|
if __name__ == "__main__":
|
|
print("e066.py: " + str(euler_066()))
|
|
assert(euler_066() == 661)
|