From 7cca53ede2ee4900fce30334fb561a24d1b5d2de Mon Sep 17 00:00:00 2001
From: Tom Smeding <hallo@tomsmeding.nl>
Date: Mon, 24 Oct 2016 14:22:24 +0200
Subject: ACTUALLY fix strongPseudoPrime2

---
 .gitignore |  1 +
 bigint.cpp | 10 ++++++++++
 bigint.h   |  1 +
 primes.cpp | 19 ++++++-------------
 primes.h   |  1 +
 5 files changed, 19 insertions(+), 13 deletions(-)

diff --git a/.gitignore b/.gitignore
index cdd3cd9..38790aa 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 *.dSYM
 cryptolib.a
 *.txt
+.DS_Store
diff --git a/bigint.cpp b/bigint.cpp
index 0d2a72b..83d00c1 100644
--- a/bigint.cpp
+++ b/bigint.cpp
@@ -542,6 +542,16 @@ bool Bigint::odd() const {
 	return !even();
 }
 
+int Bigint::trailingZeros() const {
+	if(digits.size()==0)throw domain_error("trailingZeros called on zero bigint");
+	int res=0;
+	for(digit_t dig : digits){
+		if(dig==0)res+=digit_bits;
+		else return res+__builtin_ctz(dig);
+	}
+	assert(false); //bigint is not consistent: zero while digits is not empty!
+}
+
 //Produces a string with the bytes of the mantissa in little-endian order.
 string Bigint::serialiseMantissa() const {
 	string s;
diff --git a/bigint.h b/bigint.h
index c66dcf6..0378493 100644
--- a/bigint.h
+++ b/bigint.h
@@ -100,6 +100,7 @@ public:
 	slongdigit_t lowdigits() const;
 	bool even() const;
 	bool odd() const;
+	int trailingZeros() const; //throws domain_error if called on zero value
 
 	std::string serialiseMantissa() const; //stores everything but the sign
 	void deserialiseMantissa(const std::string&); //restores non-negative number; can throw invalid_argument
diff --git a/primes.cpp b/primes.cpp
index 73e2b4c..77d547f 100644
--- a/primes.cpp
+++ b/primes.cpp
@@ -109,22 +109,15 @@ bool strongPseudoPrime2(const Bigint &n){
 	nm1-=1;
 	Bigint d(nm1);
 
-	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;
-	}
+	int zerobits=d.trailingZeros();
+	d>>=zerobits-1;
+	assert(d.even());
+	d>>=1;
+	assert(d.odd());
 	Bigint bi(expmod(Bigint::two,d,n));
 	if(bi==1||bi==nm1)return true;
 	for(int i=1;i<zerobits;i++){
-		d<<=1;
+		bi=(bi*bi).divmod(n).second;
 		if(bi==nm1)return true;
 	}
 
diff --git a/primes.h b/primes.h
index 25d0226..7d3edac 100644
--- a/primes.h
+++ b/primes.h
@@ -3,6 +3,7 @@
 #include <vector>
 #include <utility>
 #include "bigint.h"
+#include "rng.h"
 
 extern std::vector<int> smallprimes;
 
-- 
cgit v1.2.3-54-g00ecf