OFFSET

1,1

COMMENTS

Solving for z gives z = (x*y) / (x+y), so x*y == 0 (mod x+y).

All known terms are from A025487:

a(1) = 2 = 2;

a(2) = 12 = 2^2 * 3;

a(3) = 60 = 2^2 * 3 * 5;

a(4) = 840 = 2^3 * 3 * 5 * 7;

a(5) = 9240 = 2^3 * 3 * 5 * 7 * 11.

If a solution to the equation 1/z = 1/x + 1/y is found such that gcd(x,y,z) is a square, then x+y, x*y*z, and (x-y)^2 + (2*z)^2 are also squares.

For all solutions, x^2 + y^2 + z^2 is a square.

The sequence is indeed a subsequence of A025487, and likely of A126098 as well. - Max Alekseyev, Mar 01 2023

a(n) < 5*10^(n-1). - Max Alekseyev, Mar 01 2023

LINKS

Max Alekseyev, Table of n, a(n) for n = 1..30

Project Euler, Diophantine reciprocals I, Problem 108.

EXAMPLE

For n=1, we have the following, where r = (x*y) mod (x+y). (In the last four columns, each number marked by an asterisk is a square.)

.

r z x y x*y x+y x*y*z x^2+y^2+z^2

- - - - --- --- ----- -----------

0 1 2 2 4* 4* 4* 9* (solution)

2 1 2 4 8 6 8 21

4 1 2 6 12 8 12 41

6 1 2 8 16* 10 16* 69

3 1 3 3 9 6 9* 19

0 2 3 6 18* 9* 36* 49* (solution)

3 2 3 9 27 12 54 94

0 2 4 4 16* 8 32 36* (solution)

8 2 4 8 32 12 64* 84

5 2 5 5 25* 10 50 54

0 3 6 6 36* 12 108 81* (solution)

7 3 7 7 49* 14 147 107

0 4 8 8 64* 16* 256* 144* (solution)

9 4 9 9 81* 18 324* 178

.

z = 2 has the largest number of solutions, so a(1) = 2.

The number of solutions for the resulting z cannot exceed A018892(z).

PROG

(Python)

def a(n):

D = {}

for x in range(2, 10**n):

for y in range(x, 10**n):

z, r = divmod(x * y, x + y)

if r == 0:

if z in D:

D[z] += 1

else:

D[z] = 1

best_Z = 0

for z in D:

if D[z] > best_Z:

best_Z = D[z]

best_z = z

return best_z

CROSSREFS

KEYWORD

nonn

AUTHOR

Darío Clavijo, Apr 06 2022

EXTENSIONS

a(6) from Chai Wah Wu, Apr 10 2022

a(7)-a(22) from Max Alekseyev, Mar 01 2023

STATUS

approved