The OEIS mourns the passing of Jim Simons and is grateful to the Simons Foundation for its support of research in many branches of science, including the OEIS.
login
The OEIS is supported by the many generous donors to the OEIS Foundation.

 

Logo
Hints
(Greetings from The On-Line Encyclopedia of Integer Sequences!)
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. 0
0, 1, 4, 12, 40, 128, 402, 1278, 4040, 12776, 40417, 127803, 404136, 1277995, 4041401, 12779996, 40413886, 127799963 (list; graph; refs; listen; history; text; internal format)
OFFSET
1,3
LINKS
Mathematics Stack Exchange, Diophantine equation of three variables
FORMULA
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.
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...
EXAMPLE
For n = 3, the solutions are (20, 15, 12), (156, 65, 60), (136, 255, 120), (600, 175, 168). So a(3) = 4.
MATHEMATICA
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];
PROG
(Python)
from math import gcd
def fourth_root(n):
u, s = n, n + 1
while u < s:
s = u
t = 3 * s + n // (s ** 3)
u = t // 4
return s
def a(n):
N = 10 ** n
L = fourth_root(N) * 3
cnt = 0
for a in range(1, L + 1):
a2 = a * a
for b in range(1, a):
if (a + b) % 2 == 1 and gcd(a, b) == 1:
b2 = b * b
v = (4 * a * b + a2) * a2 - b2 * b2
if v > N:
continue
cnt += 1
return cnt
(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
CROSSREFS
Sequence in context: A058353 A104525 A126986 * A090576 A152174 A087206
KEYWORD
nonn,more
AUTHOR
Asif Ahmed, Feb 25 2021
STATUS
approved

Lookup | Welcome | Wiki | Register | Music | Plot 2 | Demos | Index | Browse | More | WebCam
Contribute new seq. or comment | Format | Style Sheet | Transforms | Superseeker | Recents
The OEIS Community | Maintained by The OEIS Foundation Inc.

License Agreements, Terms of Use, Privacy Policy. .

Last modified May 15 18:29 EDT 2024. Contains 372549 sequences. (Running on oeis4.)