Code:
/*
This algorithm is generic and does not exploit that q / p < 2
Plus it uses a single A and not many A
*/
A=9+16*a;//choose A with many small factors
if(M % 4 ==1){
M=3*M;
}
if((M % 8 == 3){
N=M;
}else{
N=5*M;
}
while(1){
if([1/4*(sqrt(N+1)+2)] != (int)[1/4*(sqrt(N+1)+2)]){
x=(int)[1/4*(sqrt(N+1)+2)];
}else{
x=[1/4*(sqrt(N+1)+2)]-1;
}
P=4*x+2-sqrt[4*(2*x+1)^2-N];
Q=N/P;
if (P is integer && (P % M) !=0 && (Q % M) !=0){
breack;
}
N=N*A
}
p=GCD(P,M);