cp's OEIS Frontend

A385713 Numbers w + x + y + z where (w,x,y,z) are nonnegative solutions to the Diophantine equation w^2 + x^2 + y^2 + z^2 = w*x*y*z.

%I A385713 #24 Jul 31 2025 16:56:54
%S A385713 0,8,12,32,108,292,392,1452,3392,3712,5408,11808,20172,34872,40332,
%T A385713 50572,75272,162212,280908,480512,518252,595052,700352,1048352,
%U A385713 1639308,3912492,4599172,4911172,5725732,6049292,6508932,7132608,8279808,9739812
%N A385713 Numbers w + x + y + z where (w,x,y,z) are nonnegative solutions to the Diophantine equation w^2 + x^2 + y^2 + z^2 = w*x*y*z.
%C A385713 The trivial sum with x=y=z=w=0 is included. The solutions are generated by the Vieta jumping discussed in the Numberphile video. The algorithm finds all the quadruplet solutions within the radius of a 4D sphere. Then the quadruplets are sorted based on their sum in increasing order.
%H A385713 Bence Bernáth, <a href="/A385713/b385713.txt">Table of n, a(n) for n = 1..10000</a>
%H A385713 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 A385713 8 is a term because (w,x,y,z) = (2,2,2,2) is the first nontrivial solution.
%e A385713 12 is a term because (2,2,2,6) is a solution.
%o A385713 (Python)
%o A385713 import math
%o A385713 def is_perfect_square(n):
%o A385713     return (math.isqrt(n)) ** 2 == n
%o A385713 def generate_all_solutions(radius):
%o A385713     solutions = set()
%o A385713     visited = set()
%o A385713     # Initial known solution
%o A385713     seed = (2, 2, 2, 2)
%o A385713     queue = [seed]
%o A385713     solutions.add(seed)
%o A385713     while queue:
%o A385713         quad = queue.pop(0)
%o A385713         for ii in range(4):
%o A385713             # Cyclically rotate variables: solve for the ii-th one
%o A385713             rotated = list(quad[ii:] + quad[:ii])
%o A385713             x, y, z, w = rotated
%o A385713             # Use Vieta's formula: x^2 - (yzw)x + (y² + z² + w²) = 0
%o A385713             product = y * z * w
%o A385713             sumsq = y**2 + z**2 + w**2
%o A385713             D = product**2 - 4 * sumsq
%o A385713             if D < 0 or not is_perfect_square(D):
%o A385713                 continue
%o A385713             sqrt_D = (math.isqrt(D))
%o A385713             for sign in [+1, -1]:
%o A385713                 x_new = (product + sign * sqrt_D)
%o A385713                 if x_new % 2 != 0:
%o A385713                     continue
%o A385713                 x_new //= 2
%o A385713                 if not (0 < x_new < radius):
%o A385713                     continue
%o A385713                 new_quad = [x_new, y, z, w]
%o A385713                 if any(val >= radius for val in new_quad):
%o A385713                     continue
%o A385713                 new_quad_sorted = tuple(sorted(new_quad))
%o A385713                 if new_quad_sorted not in visited:
%o A385713                     solutions.add(new_quad_sorted)
%o A385713                     queue.append(tuple(new_quad))
%o A385713                     visited.add(new_quad_sorted)
%o A385713     # Return solutions sorted by sum of elements
%o A385713     return sorted([list(sol) for sol in solutions], key=lambda tup: sum(tup))
%o A385713 print([sum(sol) for sol in generate_all_solutions(10000000)])
%Y A385713 Cf. A061292, A380462.
%K A385713 nonn
%O A385713 1,2
%A A385713 _Bence Bernáth_, Jul 07 2025