> The most difficult is to prove P is prime.

.. which takes 10^3 times longer :-)

The Cornacchia routine

{corn(p)=local(a,b,l,r);a=p;b=p-lift(sqrt(Mod(-1,p)));l=sqrtint(p);

while(b>l,r=a%b;a=b;b=r);a=sqrtint(p-b^2);[p,max(a,b),min(a,b)]}

takes Pari-GP about 2 seconds per 4*k+1 just-titanic PrP

(not fast enough for you, bien sur, but OK by me:-)

Thanks for all your guidance!

David