? f(p)={a=Mod(2^((p+1)/2)-1,2^p-1);b=a;for(i=2,p,a=sqr(a));b==a}; ? forprime(p=5,4423,if(f(p),print(p))) 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 ? ## *** last result: cpu time 3,185 ms, real time 3,186 ms. ?