121 lines
3.1 KiB
Python
121 lines
3.1 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_old(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 get_minimal_solution(d, d1=0, n1=1, d2=1, n2=1):
|
||
"""
|
||
This is an alternative to the alternate fraction approach explained here
|
||
[1]. Instead of calculating the continued fraction we search the
|
||
Stern-Brocot tree for fractions (x / y) that satisfy the pell equation.
|
||
|
||
If the solution does not satisfy the equation we can searh the respective
|
||
branch by calculating the mediant between the current fraction and one of
|
||
the previous fractions.
|
||
|
||
[1] https://en.wikipedia.org/wiki/Pell's_equation#Fundamental_solution_via_continued_fractions
|
||
[2] https://en.wikipedia.org/wiki/Stern–Brocot_tree
|
||
"""
|
||
|
||
x = n1 + n2
|
||
y = d1 + d2
|
||
p = x * x - d * y * y
|
||
if p == 1:
|
||
return (x, y)
|
||
if p < 1:
|
||
return get_minimal_solution(d, d1, n1, y, x)
|
||
else:
|
||
return get_minimal_solution(d, y, x, d2, n2)
|
||
|
||
|
||
def euler_066():
|
||
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, d_max = x, 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)
|