OFFSET
1,1
COMMENTS
EXAMPLE
a(1) = 16 = 2*(3+5).
16 is in the sequence since it is twice the sum of twin primes 3 and 5, but cannot be expressed as the sum of 2 distinct twin pairs.
36 is not in the sequence because although it is the sum of twin primes 17 and 19, it can also be expressed as the sum of pairs (5, 7) and (11, 13).
MAPLE
with(SignalProcessing): # requires at least Maple 17
N:= 10^6; # to check primes up to N
Primes:= select(isprime, {seq(2*i+1, i=1..N)}):
Twins:= Primes intersect map(t-> t-2, Primes):
nT:= nops(Twins);
T:= Array(1..(Twins[nT]+1)/2, datatype=float[8]);
for i from 1 to nT do T[(Twins[i]+1)/2]:= 1 od:
TTwins:= Convolution(T, T);
map(t -> 4*(t+1), select(n -> round(TTwins[n])=1, [$1..(nT+1)/2])); # Robert Israel, Jun 15 2014
PROG
(PARI) isok(isum1, vsum2) = {for (k=1, #vsum2, ksum2 = vsum2[k]; if (ksum2 > one, break; ); if (isum1 - ksum2 != ksum2, if (vecsearch(vsum2, isum1 - ksum2), return (0)); ); ); return (1); }
lista() = {v = readvec("b014574.txt"); vsum1 = 4*v; vsum2 = 2*v; maxs2 = vecmax(vsum2); for (i=1, #v, isum1 = vsum1[i]; if (isum1 < maxs2, if (isok(isum1, vsum2), print1(isum1, ", ")); ); ); } \\ Michel Marcus, Jun 15 2014
(PARI) l1=l2=List(); a=select(p->isprime(p+2), primes(1000)); for(i=1, #a-1, if(i<#a/4, listput(l1, 4*a[i]+4)); for(j=i+1, #a, listput(l2, 2*(a[i]+a[j])+4))); print(setminus(Set(l1), Set(l2))) \\ Lear Young, Jun 15 2014
CROSSREFS
KEYWORD
nonn
AUTHOR
Lear Young, Jun 14 2014
STATUS
approved