cp's OEIS Frontend

A341990 a(n) is the number of primitive solutions to the inverse Pythagorean equation 1/x^2 + 1/y^2 = 1/z^2 such that x <= y and x + y + z <= 10^n.

%I A341990 #39 Dec 22 2024 21:33:27
%S A341990 0,1,4,12,40,128,402,1278,4040,12776,40417,127803,404136,1277995,
%T A341990 4041401,12779996,40413886,127799963
%N A341990 a(n) is the number of primitive solutions to the inverse Pythagorean equation 1/x^2 + 1/y^2 = 1/z^2 such that x <= y and x + y + z <= 10^n.
%H A341990 Mathematics Stack Exchange, <a href="https://math.stackexchange.com/questions/2688808/diophantine-equation-of-three-variables">Diophantine equation of three variables</a>
%F A341990 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.
%F A341990 Limit_{n -> oo} a(n)/10^(n/2) = (4/Pi^2)*Integral_{x=sqrt(2)-1..1} 1/sqrt(1+4*x-x^4) ~ 0.1277999513464289283641211182341081978805...
%e A341990 For n = 3, the solutions are (20, 15, 12), (156, 65, 60), (136, 255, 120), (600, 175, 168). So a(3) = 4.
%t A341990 a[n_] := Module[{m, l, a, b, s, a2, b2, x, y, z, cnt}, m = 10^n; s = 0; l = 3 * Floor[m^0.25]; cnt = 0; For[a = 1, a <= l, a++, a2 = a * a; For[b = 1, b < a, b++, If[GCD[a, b] == 1 && Mod[a + b, 2] == 1, b2 = b * b; x = 2 * a * b * (a2 + b2); y = a2 * a2 - b2 * b2; z = 2 * a * b * (a2 - b2); If[x + y + z > m, Continue[]]; cnt += 1];]]; cnt];
%o A341990 (Python)
%o A341990 from math import gcd
%o A341990 def fourth_root(n):
%o A341990     u, s = n, n + 1
%o A341990     while u < s:
%o A341990         s = u
%o A341990         t = 3 * s + n // (s ** 3)
%o A341990         u = t // 4
%o A341990     return s
%o A341990 def a(n):
%o A341990     N = 10 ** n
%o A341990     L = fourth_root(N) * 3
%o A341990     cnt = 0
%o A341990     for a in range(1, L + 1):
%o A341990        a2 = a * a
%o A341990        for b in range(1, a):
%o A341990            if (a + b) % 2 == 1 and gcd(a, b) == 1:
%o A341990                b2 = b * b
%o A341990                v = (4 * a * b + a2) * a2 - b2 * b2
%o A341990                if v > N:
%o A341990                    continue
%o A341990                cnt += 1
%o A341990     return cnt
%o A341990 (PARI) a(n) = {my(lim = 3*sqrtnint(10^n, 4), nb = 0); for (x=1, lim, for (y=1, x, if (((x+y) % 2) && (gcd(x,y) == 1), if (2*x*y*(x^2 + y^2) + x^4 - y^4 + 2*x*y*(x^2 - y^2) <= 10^n, nb++);););); nb;} \\ _Michel Marcus_, Mar 25 2021
%K A341990 nonn,more
%O A341990 1,3
%A A341990 _Asif Ahmed_, Feb 25 2021