For Mersenne:
t1(q)={w=(2^q-1);S0=4;S=S0;for(i=1,q-1,S=Mod(S^2-2,w));s1=lift(Mod(S-5-9,w));s2=lift(Mod(S-5+9,w));if(s2==2,print1(q","))}
forprime(x=3,5000,t1(x))
Code:
3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281,3217,4253,4423,
For Wagstaff:
t2(q)={w=(2^q+1)/3;S0=4;S=S0;for(i=1,q-1,S=Mod(S^2-2,w));s1=lift(Mod(S-5-9,w));s2=lift(Mod(S-5+9,w));if((s1==0)||(s2==0),print1(q","))}
forprime(x=3,5000,t2(x))
Code:
3,5,7,11,13,17,19,23,31,43,61,79,101,127,167,191,199,313,347,701,1709,2617,3539,