A peaceful and pleasant night for you,

I have successfully implemented the following algorithm:

1. precalculate primes and n with help of f(n)=2n²-1

2. make a linear substitution with n=2*Mp*k+1 so that g(k)=8*Mp*k(Mp*k+1)+1

and calculate for the precalculated primes p the k=(n-1)*(Mp⁻¹)*(2⁻¹) mod p

3. As g(k)=1 mod Mp, I checked if the product of the sieving factors of g(k)=1 mod Mp

4. If the last condition is true, I calculate g(k) and divide by the sieving factors

and finally I checked if the remaining cofactor is a divisor of 2^Mp-1

The program seems to be right but a little bit slow ...

What a wounderful piece of work (600 lines in c), I did not parallel it.