Add (non-working yet) Lucas test

1 files changed, 69 insertions, 5 deletions

diff --git a/primes.cpp b/primes.cpp
index 135e7bb..c9d25c0 100644
--- a/primes.cpp
+++ b/primes.cpp
@@ -125,9 +125,10 @@ bool strongPseudoPrime2(const Bigint &n){
 	return false;
 }
 
-bool bailliePSW(const Bigint &n){
-	//https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test#The_test
-	if(!strongPseudoPrime2(n))return false;
+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.even())return false;
 	int D;
 	if     ((D=  5),jacobiSymbol(Bigint(D),n)==-1);
 	else if((D= -7),jacobiSymbol(Bigint(D),n)==-1);
@@ -143,6 +144,69 @@ bool bailliePSW(const Bigint &n){
 		}
 	}
 	//now we have a D
-	int P=1,Q=(1-D)/4;
-	return strongLucasPrime(n,D,P,Q);
+	// 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
+
+	//now begin the actual sequence algorithm
+#if 1
+	int s=0;
+	Bigint d(n);
+	d-=1;
+	while(d.even()){
+		d>>=1;
+		s++;
+	}
+#else
+	Bigint d(n);
+	d+=1;
+#endif
+	// cerr<<"d="<<d<<endl;
+	vector<bool> dbits=d.bits();
+	assert(dbits.size()>0);
+	assert(dbits[dbits.size()-1]==true);
+	Bigint U(1),V(P);
+	Bigint Qk(iQ);
+	for(int i=dbits.size()-2;i>=0;i--){
+		U*=V;
+		V*=V;
+		V-=Qk<<1;
+		Qk*=Qk;
+		if(dbits[i]){
+			Bigint unext(U*P);
+			unext+=V;
+			if(unext.odd())unext+=n;
+			assert(unext.even());
+			unext>>=1;
+			Bigint vnext(U*D);
+			vnext+=V*P;
+			if(vnext.odd())vnext+=n;
+			assert(vnext.even());
+			vnext>>=1;
+			U=unext.divmod(n).second;
+			V=vnext.divmod(n).second;
+			Qk=(Qk*iQ).divmod(n).second;
+		} else {
+			U=U.divmod(n).second;
+			V=V.divmod(n).second;
+			Qk=Qk.divmod(n).second;
+		}
+	}
+	if(U==0)return true;
+	if(V==0)return true; //r=0 check
+#if 1
+	for(int r=1;r<s;r++){
+		V*=V;
+		V-=Qk<<1;
+		V=V.divmod(n).second;
+		Qk=(Qk*Qk).divmod(n).second;
+		if(V==0)return true;
+	}
+#endif
+	return false;
+}
+
+bool bailliePSW(const Bigint &n){
+	//https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test#The_test
+	return strongPseudoPrime2(n)&&strongLucasPrime(n);
 }