N:= 100000: # for terms <= N
P:= select(isprime, [seq(i, i=3..2*N)]):
f:= proc(n) local m, Q, q;
m:= ListTools:-BinaryPlace(P, 2*n);
Q:= convert(P[1..m], set);
Q:= Q intersect map(t -> 2*n-t, Q);
add(2*n mod q, q = Q);
end proc:
select(t -> A[t] mod (2*t) = 0, [$1..N]);
|