revdigs:= proc(n) local L, i; L:= convert(n, base, 10); add(L[-i]*10^(i-1), i=1..nops(L)) end proc:
isemirp:= proc(p) local r;
if not isprime(p) then return false fi;
r:= revdigs(p);
r <> p and isprime(r)
end proc:
filter:= proc(p) local q, t, flag;
if not isprime(p) then return false fi;
q:= revdigs(p);
if q=p or not isprime(q) then return false fi;
flag:= false;
for t in select(`<`, numtheory:-divisors(p+1), floor(sqrt(p+1))) do
if isemirp(t-1) and isemirp((p+1)/t-1) then flag:= true; break fi
od;
if not flag then return false fi;
for t in select(`<`, numtheory:-divisors(q+1), floor(sqrt(q+1))) do
if isemirp(t-1) and isemirp((q+1)/t-1) then return true fi
od;
false
end proc:
p:= 2: R:= NULL: count:= 0:
while count < 40 do
p:= nextprime(p);
d:= ilog10(p);
p1:= floor(p/10^d);
if p1=2 then p:= nextprime(3*10^d)
elif member(p1, {4, 5, 6}) then p:= nextprime(7*10^d)
elif p1=8 then p:= nextprime(9*10^d)
fi;
if filter(p) then R:= R, p; count:= count+1 fi;
od:
R;
|