%I
%S 13,127,457,1429,1553,1621,2273,2341,2837,4129,4231,4561,4813,5119,
%T 5519,5531,6121,6451,6547,8161,8167,8219,8237,8783,8819,8831,8941,
%U 9511,10267,10559,11299,11383,12809,13183,15091,15569,16573,17569,17659,18133
%N Primes p such that p plus or minus the sum of its digits squared yields a prime in both cases.
%H Robert Israel, <a href="/A179550/b179550.txt">Table of n, a(n) for n = 1..10000</a>
%e a(5)=1553 since 1553+(1^2+5^2+5^2+3^2)=1553+60=1613 is a prime AND 1553-(1^2+5^2+5^2+3^2)=1553-60=1493 is a prime again.
%p filter:= proc(p) local t,r;
%p if not isprime(p) then return false fi;
%p r:= add(t^2, t=convert(p,base,10));
%p isprime(p+r) and isprime(p-r);
%p end proc:
%p select(filter, [seq(i,i=3..20000,2)]); # _Robert Israel_, Mar 30 2021
%t Select[Prime[Range[2100]],AllTrue[#+{Total[IntegerDigits[#]^2],-Total[ IntegerDigits[ #]^2]},PrimeQ]&] (* _Harvey P. Dale_, Aug 07 2021 *)
%o (PARI) sumdd(n) = {digs = digits(n, 10); return (sum(i=1, #digs, digs[i]^2));}
%o lista(nn) = {forprime(p=2, nn, s = sumdd(p); if (isprime(p+s) && isprime(p-s), print1(p, ", ")););} \\ _Michel Marcus_, Jul 25 2013
%o (Python)
%o from sympy import isprime, primerange
%o def sumdd(n): return sum(int(d)**2 for d in str(n))
%o def list(nn):
%o for p in primerange(2, nn+1):
%o s = sumdd(p)
%o if isprime(p-s) and isprime(p+s): print(p, end=", ")
%o list(18133) # _Michael S. Branicky_, Mar 30 2021 after _Michel Marcus_
%Y Cf. A076162, A076163, A179549.
%K nonn,base
%O 1,1
%A _Carmine Suriano_, Jul 19 2010
|