from math import atan2, cos, fmod, gcd, isqrt, pi, sin, sqrt


TARGET = 1_000_000
HALF_PI = pi / 2.0


def is_square(n):
    root = isqrt(n)
    return root if root * root == n else -1


def normalize_angle(angle):
    angle = fmod(angle, HALF_PI)
    if angle < 0.0:
        angle += HALF_PI
    if angle < 1e-18 or HALF_PI - angle < 1e-18:
        return 0.0
    return angle


def minimum_square(a, b, c):
    a, b, c = sorted((a, b, c))
    x = (a * a + c * c - b * b) / (2.0 * c)
    y2 = a * a - x * x
    if -1e-12 < y2 < 0.0:
        y2 = 0.0
    y = sqrt(y2)
    points = ((0.0, 0.0), (float(c), 0.0), (x, y))

    def evaluate(angle):
        co = cos(angle)
        si = sin(angle)
        u = [px * co + py * si for px, py in points]
        v = [-px * si + py * co for px, py in points]
        u_min = min(range(3), key=u.__getitem__)
        u_max = max(range(3), key=u.__getitem__)
        v_min = min(range(3), key=v.__getitem__)
        v_max = max(range(3), key=v.__getitem__)
        return (
            max(u[u_max] - u[u_min], v[v_max] - v[v_min]),
            u_min,
            u_max,
            v_min,
            v_max,
        )

    angles = [0.0]
    for i in range(3):
        xi, yi = points[i]
        for j in range(i + 1, 3):
            xj, yj = points[j]
            edge = atan2(yj - yi, xj - xi)
            angles.append(normalize_angle(edge))
            angles.append(normalize_angle(edge + HALF_PI))

    angles.sort()
    boundaries = []
    for angle in angles:
        if not boundaries or abs(angle - boundaries[-1]) > 1e-15:
            boundaries.append(angle)

    candidates = list(boundaries)
    endpoints = boundaries + [HALF_PI]
    for idx in range(len(endpoints) - 1):
        lo = endpoints[idx]
        hi = endpoints[idx + 1]
        if hi - lo < 1e-15:
            continue

        _, u_min, u_max, v_min, v_max = evaluate((lo + hi) / 2.0)
        ax = points[u_max][0] - points[u_min][0]
        ay = points[u_max][1] - points[u_min][1]
        bx = points[v_max][0] - points[v_min][0]
        by = points[v_max][1] - points[v_min][1]
        left = ax - by
        right = ay + bx
        if abs(right) < 1e-18:
            continue

        root = atan2(-left, right)
        for angle in (root, root + pi):
            current = normalize_angle(angle)
            if lo + 1e-15 < current < hi - 1e-15:
                candidates.append(current)

    return min(evaluate(angle)[0] for angle in candidates)


def triangle_key(a, b, c):
    a, b, c = sorted((a, b, c))
    return (a << 42) | (b << 21) | c


def pythagorean_offsets(limit):
    by_side = [[] for _ in range(limit + 1)]
    bound = isqrt(2 * limit) + 2

    for m in range(2, bound + 1):
        for n in range(1, m):
            if ((m - n) & 1) == 0 or gcd(m, n) != 1:
                continue

            a = m * m - n * n
            b = 2 * m * n
            c = m * m + n * n

            for side0, offset0 in ((a, b), (b, a)):
                side = side0
                offset = offset0
                hypotenuse = c
                while side <= limit:
                    if offset <= side:
                        by_side[side].append((offset, hypotenuse))
                    side += side0
                    offset += offset0
                    hypotenuse += c

    for side in range(limit + 1):
        if by_side[side]:
            by_side[side] = sorted(set(by_side[side]))

    return by_side


def solve(limit):
    legs = pythagorean_offsets(limit)
    seen = set()
    total = 0

    def add_triangle(a, b, c):
        nonlocal total
        key = triangle_key(a, b, c)
        if key not in seen:
            seen.add(key)
            total += a + b + c

    for side, current in enumerate(legs):
        size = len(current)
        for i in range(size):
            p, hyp_i = current[i]
            for j in range(i, size):
                q, hyp_j = current[j]
                base = p + q
                if base <= side and p * q >= side * (side - base):
                    add_triangle(base, hyp_i, hyp_j)

                u = side - p
                v = side - q
                if u > 0 and v > 0:
                    third = is_square(u * u + v * v)
                    if third >= 0 and minimum_square(hyp_i, hyp_j, third) + 1e-6 >= side:
                        add_triangle(hyp_i, hyp_j, third)

    return total


def run_checkpoints():
    assert solve(40) == 346
    assert solve(400) == 76_402
    assert solve(2_000) == 3_237_036


if __name__ == "__main__":
    run_checkpoints()
    print(solve(TARGET))