%I #72 Oct 05 2021 22:43:34
%S 6,18,54,87,108,174,174,324,324,324,492,492,492,984,984,1296,1296,
%T 1296,1440,1440,2592,2592,2592,2592,3960,3960,3960,3960,4320,4320,
%U 4320,5760,5940,5940,5940,5940,5940,5940,8640,9900,9900,9900,11880,11880,11880,11880,11880
%N a(n) is the smallest number m such that m^3 = x^3 + y^3 + z^3, x > y > z > 0, has at least n different solutions.
%C a(n) is the smallest number for which there are at least n sets of positive integers (b_i, c_i, d_i) i=1..n which satisfy the equation a(n)^3 = b_i^3 + c_i^3 + d_i^3.
%C This sequence is related to Euler's sum of powers conjecture. In particular to the case k=3, a(n) is the smallest number that has at least n different solutions to the equation.
%C The sequences of numbers whose cubes can be expressed as the sum of 3 positive cubes in at least n ways for n = 1, 2, 3, ... form a family of related sequences. This sequence is the sequence of first terms in that family of sequences.
%C The first of this family is A023042.
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Euler%27s_sum_of_powers_conjecture">Euler's sum of powers conjecture</a>
%e a(1) = 6 because 6^3 = 5^3 + 4^3 + 3^3; 6 = a(1) = A023042(1).
%e a(2) = 18 because 18^3 = 15^3 + 12^3 + 9^3 = 16^3 + 12^3 + 2^3.
%e a(3) = 54 because 54^3 = 45^3 + 36^3 + 27^3 = 48^3 + 36^3 + 6^3 = 53^3 + 19^3 + 12^3.
%o (Python)
%o import numpy as np
%o def residual(a,b,c,d, exp=3):
%o return a**exp-b**exp-c**exp-d**exp
%o def test(max_n,k=3):
%o ans=dict()
%o for a in range(max_n):
%o #print(a)
%o for b in range(int(np.ceil((a**k/3)**(1/k))),a):
%o n3=a**k-b**k
%o for c in range(int(np.ceil((n3/2)**(1/k))),b):
%o m3=n3-c**k
%o if m3<0:
%o break;
%o l=int(np.ceil((m3)**(1/k)))
%o options=[l,l-1]
%o for d in options:
%o res=residual(a,b,c,d, exp=k)
%o if res==0:
%o if a in ans.keys():
%o ans[a].append((a,b,c,d))
%o else:
%o ans[a]=[(a,b,c,d)]
%o #print("found:",(a,b,c,d))
%o break
%o else:
%o #print("tested: {0}, residual: {1}".format((a,b,c,d),res))
%o if res>0:
%o break
%o return ans
%o def serie(N):
%o result=test(N)
%o results_by_number_of_answers=[]
%o results_by_number_of_answers.append(result)
%o temp=dict()
%o for k in result.keys():
%o if len(result[k])>=2:
%o temp[k]=result[k]
%o results_by_number_of_answers.append(temp)
%o i=3
%o while len(temp)>0:
%o temp=dict()
%o for k in results_by_number_of_answers[-1].keys():
%o if len(results_by_number_of_answers[-1][k])>=i:
%o temp[k]=result[k]
%o if len(temp)>0:
%o results_by_number_of_answers.append(temp)
%o i+=1
%o return [next(iter(a)) for a in results_by_number_of_answers]
%o #Get the elements of the serie up until A_n>1000
%o A=serie(1000)
%o print(A)
%o (Python)
%o from itertools import combinations
%o from collections import Counter
%o from sympy import integer_nthroot
%o def icbrt(n): return integer_nthroot(n, 3)[0]
%o def aupto(mmax):
%o cbs = [i**3 for i in range(mmax+1)]
%o cbsset = set(cbs)
%o c = Counter(sum(c) for c in combinations(cbs, 3) if sum(c) in cbsset)
%o nmax = max(c.values())
%o return [min(icbrt(s) for s in c if c[s] >= n) for n in range(1, nmax+1)]
%o print(aupto(500)) # _Michael S. Branicky_, Sep 04 2021
%Y Cf. A023042, A025418, A346137, A316359.
%K nonn,more
%O 1,1
%A _Sebastian Magee_, Jul 30 2021
%E a(16)-a(31) from _Jinyuan Wang_, Aug 02 2021
%E More terms from _David A. Corneth_, Sep 04 2021