Add (non-working yet) Lucas test

5 files changed, 94 insertions, 12 deletions

diff --git a/bigint.cpp b/bigint.cpp
index 21c4ab2..9fdafe0 100644
--- a/bigint.cpp
+++ b/bigint.cpp
@@ -133,8 +133,9 @@ Bigint& Bigint::operator=(slongdigit_t v){
 }
 
 Bigint& Bigint::operator+=(const Bigint &o){
-	if(&o==this){
-		return *this=Bigint(*this)+=o;
+	if(&o==this){ //*this + *this
+		operator<<=(1);
+		return *this;
 	}
 	if(sign==1){
 		if(o.sign==1)add(*this,o);
@@ -148,8 +149,10 @@ Bigint& Bigint::operator+=(const Bigint &o){
 }
 
 Bigint& Bigint::operator-=(const Bigint &o){
-	if(&o==this){
-		return *this=Bigint(*this)-=o;
+	if(&o==this){ // *this - *this
+		sign=1;
+		digits.clear();
+		return *this;
 	}
 	if(sign==1){
 		if(o.sign==1)subtract(*this,o);
diff --git a/main.cpp b/main.cpp
index db846a9..43d2de7 100644
--- a/main.cpp
+++ b/main.cpp
@@ -153,6 +153,16 @@ void fermatpseudo(){
 	cout<<endl;
 }
 
+void lucaspseudo(){
+	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<<' ';
+	}
+	cout<<endl;
+}
+
 int main(int argc,char **argv){
 	(void)argc;
 	(void)argv;
@@ -162,4 +172,6 @@ int main(int argc,char **argv){
 	// testisqrt(argc,argv);
 	// randprime(Bigint(20),Bigint(42));
 	// fermatpseudo();
+	// strongLucasPrime(Bigint(5));
+	lucaspseudo();
 }
diff --git a/numalgo.cpp b/numalgo.cpp
index 9b1681c..c3bfa1b 100644
--- a/numalgo.cpp
+++ b/numalgo.cpp
@@ -8,7 +8,7 @@ Bigint gcd(Bigint a,Bigint b){
 	while(true){
 		if(a==0)return b;
 		if(b==0)return a;
-		if(a>=b)a=a.divmod(b).second;
+		if(a.compareAbs(b)==1)a=a.divmod(b).second;
 		else b=b.divmod(a).second;
 	}
 }
@@ -45,13 +45,16 @@ Bigint expmod(const Bigint &b,const Bigint &e,const Bigint &m){
 
 int jacobiSymbol(Bigint a,Bigint n){
 	//https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol
+	// cerr<<"jacobiSymbol("<<a<<','<<n<<')'<<endl;
 	assert(n>0&&n.odd());
 	int res=1;
 	while(true){
+		// cerr<<" a="<<a<<" n="<<n<<endl;
 		a=a.divmod(n).second;
 		if(a<0)a+=n;
 		int nmod8=n.divmod(Bigint(8)).second.lowdigits();
 		int two_over_n=nmod8==1||nmod8==7?1:-1;
+		if(a==0)return 0;
 		while(a.even()){
 			res*=two_over_n;
 			a>>=1;
@@ -84,7 +87,7 @@ Bigint isqrt(const Bigint &n){
 		if(xnext==x)break;
 		x=xnext;
 	}
-	cerr<<iter<<endl;
+	// cerr<<iter<<endl;
 	assert(iter<maxiter);
 	switch((x*x).compare(n)){
 		case -1:{
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);
 }
diff --git a/primes.h b/primes.h
index 744fa3e..eb817b1 100644
--- a/primes.h
+++ b/primes.h
@@ -19,7 +19,7 @@ Bigint randprime(const Bigint &low,const Bigint &high);
 bool strongPseudoPrime2(const Bigint &n);
 
 //checks Lucas [pseudo- or actual] primality
-bool strongLucasPrime(const Bigint &n,int D,int P,int Q);
+bool strongLucasPrime(const Bigint &n);
 
 //checks primality
 bool bailliePSW(const Bigint &n);