Philmoore's algorithm with GP/Pari.
 

Here is the GP code for Philmoore's algorithm (for Mersenne primes), tested with gp v 2.1.1 .
It's quite slow! (I'm not using GMP kernel).
[code:1]
tony(q) =
{
p=2^q-1;
b=3^(2^(q-2)) % p;
if(2*b>p,b=p-b);
a=p;
s=floor(sqrt(p));
while(b>s,t=a % b; a=b; b=t);
x=b;
y=floor(sqrt((p-x^2)/3));
print("Mq= ",p);
print("x= ",x);
print("y= ",y);
print("x^2+3y^2= ",x^2+3*y^2);
}
[/code:1]

How to use:
? tony(19)
Mq= 524287
x= 298
y= 381
x^2+3y^2= 524287