cp's OEIS Frontend

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.

%I A380462 #37 Jul 01 2025 10:53:55
%S A380462 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,
%T A380462 86,94,101,111,118,132,141,150,159,168,180,192,203,220,235,247,261,
%U A380462 275,289,302,324,340,361,376,394,413,433,454,472,498,518,536,560,584,606,632,658,684
%N 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.
%C A380462 The trivial x=y=z=w=0 solution is not included. The solutions are generated by the Vieta jumping discussed in the Numberphile video.
%H A380462 Bence Bernáth, <a href="/A380462/b380462.txt">Table of n, a(n) for n = 1..650</a>
%H A380462 Michael Magee and Brady Haran, <a href="https://youtube.com/watch?v=a7BVL1MOCl4">A simple equation that behaves weirdly</a>, Numberphile, Youtube video, 2025.
%e A380462 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.
%o A380462 (Python)
%o A380462 import math
%o A380462 from collections import deque
%o A380462 def is_perfect_square(n):
%o A380462     return (math.isqrt(n)) ** 2 == n
%o A380462 def generate_all_solutions(up_to_bound):
%o A380462     solutions = set()
%o A380462     visited = set()
%o A380462     seed = (2, 2, 2, 2)
%o A380462     queue = deque([seed])
%o A380462     solutions.add(seed)
%o A380462     while queue:
%o A380462         quad = queue.popleft()
%o A380462         for ii in range(4):
%o A380462             rotated = list(quad[ii:] + quad[:ii])
%o A380462             x, y, z, w = rotated
%o A380462             product = y * z * w
%o A380462             sumsq = y**2 + z**2 + w**2
%o A380462             D = product**2 - 4 * sumsq
%o A380462             if D < 0 or not is_perfect_square(D):
%o A380462                 continue
%o A380462             sqrt_D = (math.isqrt(D))
%o A380462             for sign in [+1, -1]:
%o A380462                 x_new = (product + sign * sqrt_D)
%o A380462                 if x_new % 2 != 0:
%o A380462                     continue
%o A380462                 x_new //= 2
%o A380462                 if not (0 < x_new < up_to_bound):
%o A380462                     continue
%o A380462                 new_quad = [x_new, y, z, w]
%o A380462                 if any(val >= up_to_bound for val in new_quad):
%o A380462                     continue
%o A380462                 new_quad_sorted = tuple(sorted(new_quad))
%o A380462                 if new_quad_sorted not in visited:
%o A380462                     solutions.add(new_quad_sorted)
%o A380462                     queue.append(tuple(new_quad))
%o A380462                     visited.add(new_quad_sorted)
%o A380462     return sorted(solutions)
%o A380462 # Compute all solutions up to e^100
%o A380462 max_bound = math.exp(100)
%o A380462 all_solutions = generate_all_solutions(max_bound)
%o A380462 # Precompute max entry of each solution for fast thresholding
%o A380462 solution_maxes = [max(sol) for sol in all_solutions]
%o A380462 # Generate the sequence a(n) for n = 1 to 100
%o A380462 a_n = []
%o A380462 for n in range(1, 101):
%o A380462     threshold = math.exp(n)
%o A380462     count = sum(1 for m in solution_maxes if m < threshold)
%o A380462     a_n.append(count)
%o A380462 print(a_n)
%Y A380462 Cf. A061292.
%K A380462 nonn
%O A380462 1,2
%A A380462 _Bence Bernáth_, Jun 22 2025