login
The OEIS Foundation is supported by donations from users of the OEIS and by a grant from the Simons 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

Table of n, a(n) for n=1..18.

Math 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.

Lim_{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

Adjacent sequences:  A341987 A341988 A341989 * A341991 A341992 A341993

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 | Recent
The OEIS Community | Maintained by The OEIS Foundation Inc.

License Agreements, Terms of Use, Privacy Policy. .

Last modified July 23 14:44 EDT 2021. Contains 346259 sequences. (Running on oeis4.)