After bzip2 -d u.gp.bz2, executimg the 171MB PARI/GP script computes integers x,y from the stored sqrt(Mod(-3,M_i)) and verifies that M_i==x^2+3*y^2 in less than two minutes for the top 20 known Mersenne primes:
hermann@7950x:~/M_33-52$ time ( gp -q < u.gp )
859433: 1
1257787: 1
1398269: 1
2976221: 1
3021377: 1
6972593: 1
13466917: 1
20996011: 1
24036583: 1
25964951: 1
30402457: 1
32582657: 1
37156667: 1
42643801: 1
43112609: 1
57885161: 1
74207281: 1
77232917: 1
82589933: 1
136279841: 1
real 1m14.267s
user 1m13.777s
sys 0m0.484s
hermann@7950x:~/M_33-52$
PARI/GP script is big, but hs only few lines:
hermann@7950x:~/M_33-52$ cut -b-50 u.gp
{a=[
[1,859433,0x94e81c958670c55396eb18099f0d2f693deb63
[1,1257787,0x34d977ff67ced94b8362a29fccdfc269cb51d
[1,1398269,0x1a5bfb4f9b1077db3012c0965e9ee870e8881
[1,2976221,0x43c22f172e170669d2770027991d48c3963a7
[1,3021377,0xa68a4d9631e1e9018f8be6cea4271f9e5fccb
[1,6972593,0x1d775fde1cc545783f5d9c8fc43882618ae28
[1,13466917,0x1123b1f5d4d113e93722ab4afae6bb9695fd
[1,20996011,0x594d870b29c3650aa571ff8d9785517b6c00
[1,24036583,0x6e2411ec3042c6ab4e5b3357f50dab2dae51
[1,25964951,0xf8dc239be0db7926bf486de6a93e9c49ecdd
[1,30402457,0x15c534ca70f527658cd4ba2e378f16b042b4
[1,32582657,0x1ca57bf12727f276a4ab6e47500db92bc41f
[1,37156667,0x89682eb54f0e692948a4e47af480bd2a6be5
[1,42643801,0x156e0ee8a3dc41a50a2ae1f01f1c26c170ab
[1,43112609,0x931875a80ea54c22a14fc58a99ed0471bfad
[1,57885161,0x980ed5371fa644024f0132a6820826ea5de8
[1,74207281,0x128eaf76450bc3352e2cc367383902eaefcf
[1,77232917,0xe00f9eb438de590c2c3d9560063d81f3fc7f
[1,82589933,0xc50ae183b02da62218fcb62298eda4403a2d
[1,136279841,0x842168d36153625461ffac45c977429bd2a
];
foreach(a,t,
p=2^t[2]-1;
[M,V]=halfgcd(t[3],p);
[x,y]=[V[2],M[2,1]];
print(t[2],": ",x^2+3*y^2==p)
)}
hermann@7950x:~/M_33-52$