From 62dab77c1769c398dfb3f0a6968e35f63869e7e5 Mon Sep 17 00:00:00 2001
From: tomsmeding <tom.smeding@gmail.com>
Date: Wed, 5 Oct 2016 09:53:50 +0200
Subject: Fix Baillie-PSW!

---
 main.cpp   | 27 +++++++++++++--------------
 primes.cpp | 15 ++++++++++-----
 2 files changed, 23 insertions(+), 19 deletions(-)

diff --git a/main.cpp b/main.cpp
index 43d2de7..380fd0f 100644
--- a/main.cpp
+++ b/main.cpp
@@ -144,21 +144,19 @@ void performrsa(){
 	cout<<"msg = "<<msg2<<endl;
 }
 
-void fermatpseudo(){
-	fillsmallprimes();
-	for(int i=2;i<65000;i++){
-		if(!strongPseudoPrime2(Bigint(i)))continue;
-		if(!binary_search(smallprimes.begin(),smallprimes.end(),i))cout<<i<<' ';
-	}
-	cout<<endl;
-}
-
-void lucaspseudo(){
+void pseudolist(bool(*func)(const Bigint&)){
 	fillsmallprimes();
 	for(int i=2;i<65000;i++){
 		// cerr<<"TRYING "<<i<<endl;
-		if(!strongLucasPrime(Bigint(i)))continue;
-		if(!binary_search(smallprimes.begin(),smallprimes.end(),i))cout<<i<<' ';
+		bool actualprime=binary_search(smallprimes.begin(),smallprimes.end(),i);
+		if(!func(Bigint(i))){
+			if(actualprime){
+				cerr<<"Test misleadingly said that "<<i<<" isn't prime, while it is!"<<endl;
+				exit(1);
+			}
+			continue;
+		}
+		if(!actualprime)cout<<i<<' ';
 	}
 	cout<<endl;
 }
@@ -171,7 +169,8 @@ int main(int argc,char **argv){
 	// performrsa();
 	// testisqrt(argc,argv);
 	// randprime(Bigint(20),Bigint(42));
-	// fermatpseudo();
 	// strongLucasPrime(Bigint(5));
-	lucaspseudo();
+	// cout<<strongLucasPrime(Bigint(5777))<<endl;
+	// pseudolist(strongPseudoPrime2);
+	pseudolist(strongLucasPrime);
 }
diff --git a/primes.cpp b/primes.cpp
index c9d25c0..d259939 100644
--- a/primes.cpp
+++ b/primes.cpp
@@ -113,7 +113,9 @@ Bigint randprime(const Bigint &biglow,const Bigint &bighigh){
 
 bool strongPseudoPrime2(const Bigint &n){
 	//https://en.wikipedia.org/wiki/Strong_pseudoprime#Formal_definition
-	if(n<2||n.even())return false;
+	if(n<2)return false;
+	if(n==2)return true;
+	if(n.even())return false;
 	Bigint nm1(n);
 	nm1-=1;
 	Bigint d(nm1);
@@ -128,6 +130,8 @@ bool strongPseudoPrime2(const Bigint &n){
 bool strongLucasPrime(const Bigint &n){
 	//https://en.wikipedia.org/wiki/Lucas_pseudoprime#Implementing_a_Lucas_probable_prime_test
 	//TODO: use d and s already calculated in strongPseudoPrime2
+	if(n<2)return false;
+	if(n==2)return true;
 	if(n.even())return false;
 	int D;
 	if     ((D=  5),jacobiSymbol(Bigint(D),n)==-1);
@@ -146,13 +150,14 @@ bool strongLucasPrime(const Bigint &n){
 	//now we have a D
 	// cerr<<"D="<<D<<endl;
 	int P=1,iQ=(1-D)/4;
-	if(gcd(n,Bigint(P))!=1||gcd(n,Bigint(iQ))!=1)return false; //would be too easy to ignore
+	if(gcd(n,Bigint(P)).compareAbs(1)!=0||gcd(n,Bigint(iQ)).compareAbs(1)!=0)return false; //would be too easy to ignore
+	// cerr<<"ok"<<endl;
 
 	//now begin the actual sequence algorithm
-#if 1
+#if 0
 	int s=0;
 	Bigint d(n);
-	d-=1;
+	d+=1;
 	while(d.even()){
 		d>>=1;
 		s++;
@@ -193,8 +198,8 @@ bool strongLucasPrime(const Bigint &n){
 		}
 	}
 	if(U==0)return true;
+#if 0
 	if(V==0)return true; //r=0 check
-#if 1
 	for(int r=1;r<s;r++){
 		V*=V;
 		V-=Qk<<1;
-- 
cgit v1.2.3-70-g09d2