1 files changed, 19 insertions, 4 deletions
diff --git a/primes.cpp b/primes.cpp
index e0a7e59..73e2b4c 100644
--- a/primes.cpp
+++ b/primes.cpp
@@ -108,11 +108,26 @@ bool strongPseudoPrime2(const Bigint &n){
 	Bigint nm1(n);
 	nm1-=1;
 	Bigint d(nm1);
-	while(d.even()){ //TODO: optimise using __builtin_ctz
-		d>>=1;
-		if(expmod(Bigint::two,d,n)==nm1)return true;
+
+	int zerobits=0;
+	while(d.even()){
+		Bigint::slongdigit_t lowdig=d.lowdigits();
+		int trzeros;
+		if(lowdig==0){
+			trzeros=8*sizeof(Bigint::slongdigit_t)-1;
+		} else {
+			trzeros=__builtin_ctz(d.lowdigits());
+		}
+		zerobits+=trzeros;
+		d>>=trzeros;
+	}
+	Bigint bi(expmod(Bigint::two,d,n));
+	if(bi==1||bi==nm1)return true;
+	for(int i=1;i<zerobits;i++){
+		d<<=1;
+		if(bi==nm1)return true;
 	}
-	if(expmod(Bigint::two,d,n)==1)return true;
+
 	return false;
 }