A380462 a(n) is the number of positive solutions to the Diophantine equation w^2 + x^2 + y^2 + z^2 = w*x*y*z such that x,y,z,w < e^n.

1, 2, 2, 3, 4, 6, 6, 7, 10, 12, 16, 17, 18, 23, 25, 32, 35, 39, 43, 47, 55, 60, 64, 73, 76, 86, 94, 101, 111, 118, 132, 141, 150, 159, 168, 180, 192, 203, 220, 235, 247, 261, 275, 289, 302, 324, 340, 361, 376, 394, 413, 433, 454, 472, 498, 518, 536, 560, 584, 606, 632, 658, 684

Offset: 1

Author

Bence Bernáth, Jun 22 2025

Keywords

• nonn

Comments

The trivial x=y=z=w=0 solution is not included. The solutions are generated by the Vieta jumping discussed in the Numberphile video.

Examples

			a(3)=2 because the solutions below e^3=20.085... are [2,2,2,2] and [2,2,2,6]. One quadruplet solution is counted only once.
		

Links

• Bence Bernáth, Table of n, a(n) for n = 1..650
• Michael Magee and Brady Haran, A simple equation that behaves weirdly, Numberphile, Youtube video, 2025.

Crossrefs

Cf. A061292.

Programs

Python

import math
from collections import deque
def is_perfect_square(n):
    return (math.isqrt(n)) ** 2 == n
def generate_all_solutions(up_to_bound):
    solutions = set()
    visited = set()
    seed = (2, 2, 2, 2)
    queue = deque([seed])
    solutions.add(seed)
    while queue:
        quad = queue.popleft()
        for ii in range(4):
            rotated = list(quad[ii:] + quad[:ii])
            x, y, z, w = rotated
            product = y * z * w
            sumsq = y**2 + z**2 + w**2
            D = product**2 - 4 * sumsq
            if D < 0 or not is_perfect_square(D):
                continue
            sqrt_D = (math.isqrt(D))
            for sign in [+1, -1]:
                x_new = (product + sign * sqrt_D)
                if x_new % 2 != 0:
                    continue
                x_new //= 2
                if not (0 < x_new < up_to_bound):
                    continue
                new_quad = [x_new, y, z, w]
                if any(val >= up_to_bound for val in new_quad):
                    continue
                new_quad_sorted = tuple(sorted(new_quad))
                if new_quad_sorted not in visited:
                    solutions.add(new_quad_sorted)
                    queue.append(tuple(new_quad))
                    visited.add(new_quad_sorted)
    return sorted(solutions)
# Compute all solutions up to e^100
max_bound = math.exp(100)
all_solutions = generate_all_solutions(max_bound)
# Precompute max entry of each solution for fast thresholding
solution_maxes = [max(sol) for sol in all_solutions]
# Generate the sequence a(n) for n = 1 to 100
a_n = []
for n in range(1, 101):
    threshold = math.exp(n)
    count = sum(1 for m in solution_maxes if m < threshold)
    a_n.append(count)
print(a_n)