OFFSET
1,1
LINKS
Robert Israel, Table of n, a(n) for n = 1..564
EXAMPLE
a(1) = 6 as 6 = 2 + 2^2.
a(2) = 11 as 11 = 7 + 2^2 = 2 + 3^2.
a(3) = 56 as 56 = 47 + 3^2 = 31 + 5^2 = 7 + 7^2.
a(4) = 176 as 176 = 167 + 3^2 = 151 + 5^2 = 127 + 7^2 = 7 + 13^2.
a(5) = 188 as 188 = 179 + 3^2 = 163 + 5^2 = 139 + 7^2 = 67 + 11^2 = 19 + 13^2.
a(6) = 362 as 362 = 353 + 3^2 = 337 + 5^2 = 313 + 7^2 = 241 + 11^2 = 193 + 13^2 = 73 + 17^2.
a(9) = 1568 as 1568 = 1559 + 3^2 = 1543 + 5^2 = 1447 + 11^2 = 1399 + 13^2 = 1279 + 17^2 = 1039 + 23^2 = 727 + 29^2 = 607 + 31^2 = 199 + 37^2.
Note that a(9) > A381333(9) = 1448 as 1448 has 10 decompositions: 1448 = 1439 + 3^2 = 1423 + 5^2 = 1399 + 7^2 = 1327 + 11^2 = 1279 + 13^2 = 1087 + 19^2 = 919 + 23^2 = 607 + 29^2 = 487 + 31^2 = 79 + 37^2.
MAPLE
N:= 10^6: # for a(1)..a(k) where a(k+1) is the first term > N
P:= select(isprime, [2, seq(i, i=3..N, 2)]):
W:= Vector(1..N/2, datatype=integer[4]):
for i from 2 while P[i]^2 < N do
m:= ListTools:-BinaryPlace(P, N - P[i]^2);
J:= (P[2..m] +~ P[i]^2)/~ 2;
W[J]:= W[J] +~ 1;
od:
imax:= max[index](W):
R:= Vector(W[imax]):
R[1]:= 6: R[2]:= 11:
for i from 1 to imax do
r:= W[i];
if r > 0 and R[r] = 0 then R[r]:= 2*i fi;
od:
if member(0, R, 'i') then convert(R[1..i-1], list) else convert(R, list) fi; # Robert Israel, Feb 24 2025
PROG
(Python)
from itertools import count
from math import isqrt
from sympy import isprime, primerange
def A381334(n): return next(filter(lambda m:sum(1 for p in primerange(isqrt(m)+1) if isprime(m-p**2))==n, count(1)))
CROSSREFS
KEYWORD
nonn
AUTHOR
Chai Wah Wu, Feb 20 2025
STATUS
approved