solve_reciprocal_pythagorean_equation.sf

#!/usr/bin/ruby

# Find all the primitive solutions to the inverse Pythagorean equation:
#   1/x^2 + 1/y^2 = 1/z^2
# such that x <= y and 1 <= x,y,z <= N.

# It can be shown that all the primitive solutions are generated from the following parametric form:
#
#   x = 2*a*b*(a^2 + b^2)
#   y = a^4 - b^4
#   z = 2*a*b*(a^2 - b^2)
#
# where gcd(a, b) = 1 and a+b is odd.

# See also:
#   https://oeis.org/A341990
#   https://math.stackexchange.com/questions/2688808/diophantine-equation-of-three-variables

func S(N) {

    var solutions = []
    var limit = N.iroot(4)

    for a in (1..limit), b in (1 ..^ a) {

        is_odd(a+b) || next
        is_coprime(a,b) || next

        var x = (2*a*b * (a**2 + b**2))
        var y = (a**4 - b**4)
        var z = (2*a*b * (a**2 - b**2))

        x <= N || next
        y <= N || next
        z <= N || next

        solutions << [x,y,z]
        assert_eq(1/x**2 + 1/y**2, 1/z**2)
    }

    solutions.sort
}

var N = 1e4

say <<"EOT"

:: Primitve solutions (x,y,z) with 1 <= x,y,z <= #{N} and x <= y, to equation:

    1/x^2 + 1/y^2 = 1/z^2
EOT

for x,y,z in (S(N)) {
    say "(#{x}, #{y}, #{z})"
}

__END__

:: Primitve solutions (x,y,z) with 1 <= x,y,z <= 10000 and x <= y, to equation:

    1/x^2 + 1/y^2 = 1/z^2

(20, 15, 12)
(136, 255, 120)
(156, 65, 60)
(444, 1295, 420)
(580, 609, 420)
(600, 175, 168)
(1040, 4095, 1008)
(1484, 2385, 1260)
(1640, 369, 360)
(2020, 9999, 1980)
(3060, 6545, 2772)
(3504, 4015, 2640)
(3640, 2145, 1848)
(3660, 671, 660)
(6540, 9919, 5460)
(6984, 6305, 4680)
(7120, 3471, 3120)
(7140, 1105, 1092)