euler/python/e075.py

from math import sqrt
from lib_misc import gcd


def proper_divisors(n):
    """
    Returns the list of divisors for n excluding n.
    """
    if n < 2:
        return []
    divisors = [1, ]
    d = 2
    while d * d <= n:
        if n % d == 0:
            divisors.append(d)
        d += 1
    return divisors


def triples_brute_force(b_max):
    squares = {n * n: n for n in range(1, b_max)}
    triples = []

    for b in range(3, b_max):
        for a in range(2, b):
            c_squared = a * a + b * b
            if c_squared in squares:
                triples.append((a, b, squares[c_squared]))
    return triples


def dicksons_method(r):
    triples = []
    assert(r % 2 == 0)
    st = r * r // 2

    for d in proper_divisors(st):
        s = d
        t = st // s

        a = r + s
        b = r + t
        c = r + s + t
        L = a + b + c
        triples.append(L)
    return sorted(triples)


def dickson_method_min_max_perimeter(r):
    """ I used this method to check that the min and max
    perimeter for increasing r are also increasing strictly
    monotonic. """
    assert(r % 2 == 0)
    st = r * r // 2
    L_min = 10**8
    L_max = 0

    for d in proper_divisors(st):
        s = d
        t = st // s

        L = r + s + r + t + r + s + t
        if L > L_max:
            L_max = L
        if L < L_min:
            L_min = L
    return (L_min, L_max)


def dicksons_method_efficient(r):
    L_max = 1500000
    perimeters = []
    assert(r % 2 == 0)
    st = r * r // 2
    s_max = int(sqrt(st))
    for s in range(s_max, 0, -1):
        if st % s == 0:
            t = st // s
            L = r + s + r + t + r + s + t
            if L > L_max:
                break
            perimeters.append(L)
        else:
            pass
    return perimeters


def euler_075_dicksons():
    """
    We are using Dickson's method to get all Pythagorean triples.  The
    challenge is that this method requires finding the divisors for a
    given r. To make it more efficient we stop iterating after the
    perimeter exceeds the threshold L_max for any given r.
    In addition we found a r_max of 258000 empirically.
    This runs for 15 minutes with cpython and in under 5 minutes in pypy.
    It would have been better to use a generator function that finds native
    solutions and then extrapolate with factors k in range(2, 1500000)
    while L < L_max.

    https://en.wikipedia.org/wiki/Formulas_for_generating_Pythagorean_triples#Dickson's_method   

    XXX: Can be strongly optimized.
    """

    r_max = 258000
    perimeter_counts = {}

    for r in range(2, r_max, 2):
        if r % 10000 == 0:
            print(r)
        for perimeter in dicksons_method_efficient(r):
            try:
                perimeter_counts[perimeter] += 1
            except KeyError:
                perimeter_counts[perimeter] = 1

    one_integer_sided_right_triangle_count = 0
    for perimeter, count in perimeter_counts.items():
        if count == 1:
            one_integer_sided_right_triangle_count += 1
    return one_integer_sided_right_triangle_count


def triples_euclids_formula(m, n):
    """
    Generates only primitive triples.
    https://en.wikipedia.org/wiki/Pythagorean_triple#Generating_a_triple
    """
    assert(m > n)
    assert(n > 0)
    assert(gcd(m, n) == 1)
    assert(m % 2 == 0 or n % 2 == 0)
    a = m * m - n * n
    b = 2 * m * n
    c = m * m + n * n
    if a < b:
        return (a, b, c)
    else:
        return (b, a, c)


def euler_075():
    """
    The range for m was chosen empirically (aka by trial and error).
    I tried terminating when L > L_max the first time, but there might
    be combinations for m = m + 1 that yield lower L. I guess the right
    criteria to return would be if L > L_max and n == 1.
    This is good enough for now. I still like my other solution, but it
    is so much slower because finding the divisors is so expensive.
    """
    L = 0
    L_max = 1500000
    m = 1
    perimeter_counts = {}
    for m in range(1, 1000):
        for n in range(1, m):
            try:
                L = sum(triples_euclids_formula(m, n))
                L_k = L
                while L_k <= L_max:
                    try:
                        perimeter_counts[L_k] += 1
                    except KeyError:
                        perimeter_counts[L_k] = 1
                    L_k += L
            except AssertionError:
                pass
    one_integer_sided_right_triangle_count = 0
    for perimeter, count in perimeter_counts.items():
        if count == 1:
            one_integer_sided_right_triangle_count += 1
    return one_integer_sided_right_triangle_count


if __name__ == "__main__":
    print("e075.py: " + str(euler_075()))
    assert(euler_075() == 161667)