|
%I #18 Aug 28 2018 14:59:29
%S 2,6,5,13,9,20,30,19,52,78,18,7,43,151,79,126,88,373,183,84,177,521,
%T 263,2347,1305,392,294,1207,3946,1973,2099,185,999,518,1970,9791,4577,
%U 6111,4811,21372,10154,3210,9874,89482,49678,9918,8344,46684,65588,50136
%N a(1) = 2; for n > 1, a(n) is the least positive number not yet in the sequence such that Sum_{k=1..n} a(k) divides Sum_{k=1..n} a(k)^2.
%C Is this sequence infinite?
%C The variant of this sequence starting with 1 has only one term.
%C See A318359 for a similar sequence.
%H Chai Wah Wu, <a href="/A318358/b318358.txt">Table of n, a(n) for n = 1..250</a> (n = 1..100 from Rémy Sigrist)
%e For n = 3:
%e - (2^2 + 6^2 + 1^2) == 5 mod (2 + 6 + 1),
%e - (2^2 + 6^2 + 3^2) == 5 mod (2 + 6 + 3),
%e - (2^2 + 6^2 + 4^2) == 8 mod (2 + 6 + 4),
%e - (2^2 + 6^2 + 5^2) == 0 mod (2 + 6 + 5),
%e - hence a(3) = 5.
%o (PARI) s=0; s2=0; p=0; v=2; for (n=1, 50, print1 (v ", "); s+=v; s2+=v^2; p+=2^v; for (w=1, oo, if (!bittest(p,w) && (s2+w^2)%(s+w)==0, v=w; break)))
%o (Python)
%o import bisect
%o from sympy.solvers.diophantine import diop_quadratic
%o from sympy.abc import x,y
%o A318358_list, A318358_set, p, q = [2], {2}, 2, 2**2
%o for _ in range(100):
%o r = sorted(next(zip(*diop_quadratic(x**2+q-p*y-x*y))))
%o for a in r[bisect.bisect_right(r,0):]:
%o if a not in A318358_set:
%o A318358_list.append(a)
%o A318358_set.add(a)
%o break
%o p += a
%o q += a**2 # _Chai Wah Wu_, Aug 28 2018
%Y Cf. A318359.
%K nonn
%O 1,1
%A _Rémy Sigrist_, Aug 24 2018
|
