filter:= proc(p) local r, q;
r:= numtheory:-tau((p-1)/2^padic:-ordp(p-1, 2));
if not isprime(r) then return false fi;
r = numtheory:-tau((p+1)/2^padic:-ordp(p+1, 2))
end proc:
res:= NULL: p:= 0:
while p < 1000 do
p:= nextprime(p);
if filter(p) then
res:= res, p;
fi;
od: