|
(PARI) forprime(i=1, 10^12, di=digits(i); if(di==Vecrev(di), s=sum(j=1, #di, di[j]!); ds=digits(s); if(isprime(s) && ds==Vecrev(ds), print1(i, ", "))))
(Python)
from math import factorial
from sympy import isprime
from itertools import count, islice, product
myfact = {d:factorial(int(d)) for d in "0123456789"}
def ispal(s): return s == s[::-1]
def sofod(s): return sum(myfact[d] for d in s)
def palprimecandstrs():
yield from ["2", "3", "5", "7", "11"]
for digs in count(3, step=2):
for last in "1379":
for p in product("0123456789", repeat=digs//2-1):
left = "".join(p)
for mid in "01": # else sofod even (cf. Corneth comment)
yield last + left + mid + left[::-1] + last
def agen(): # generator of terms
for strp in palprimecandstrs():
s = sofod(strp)
if ispal(str(s)) and isprime(s):
p = int(strp)
if isprime(p):
yield p
| 