euler/python/e066.py

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)