More primes and primality testing

1 files changed, 68 insertions, 2 deletions

diff --git a/primes.cpp b/primes.cpp
index d144e5e..135e7bb 100644
--- a/primes.cpp
+++ b/primes.cpp
@@ -1,3 +1,4 @@
+#include <stdexcept>
 #include <cstring>
 #include <cmath>
 #include <cassert>
@@ -10,6 +11,7 @@ vector<int> smallprimes;
 bool smallprimes_inited=false;
 
 void fillsmallprimes(){
+	// cerr<<"Filling small primes..."<<endl;
 	smallprimes_inited=true;
 	//TODO: reserve expected amount of space in smallprimes
 	smallprimes.push_back(2);
@@ -44,6 +46,7 @@ pair<Bigint,Bigint> genprimepair(int nbits){
 }
 
 Bigint randprime(const Bigint &biglow,const Bigint &bighigh){
+	//https://en.wikipedia.org/wiki/Generating_primes#Large_primes
 	if(!smallprimes_inited)fillsmallprimes();
 
 	assert(bighigh>biglow);
@@ -53,30 +56,93 @@ Bigint randprime(const Bigint &biglow,const Bigint &bighigh){
 	if(diff<=maxrangesize){
 		low=biglow;
 		high=bighigh;
+		// cerr<<"low=biglow="<<low<<" high=bighigh="<<high<<endl;
 	} else {
 		high=low=cryptrandom_big(diff-maxrangesize);
 		high+=maxrangesize;
+		// cerr<<"low="<<low<<" high="<<high<<endl;
 	}
 	if(low.even())low+=1;
 	if(high.even())high-=1;
+	// cerr<<"low="<<low<<" high="<<high<<endl;
 	Bigint sizeb(high-low+1);
 	assert(sizeb>=0&&sizeb<=maxrangesize);
 	int nnums=(sizeb.lowdigits()+1)/2;
+	// cerr<<"nnums="<<nnums<<endl;
+
+	int nleft=nnums;
 
 	bool composite[nnums]; //low, low+2, low+4, ..., high (all odd numbers)
 	memset(composite,0,nnums*sizeof(bool));
 	int nsmallprimes=smallprimes.size();
 	for(int i=1;i<nsmallprimes;i++){
 		int pr=smallprimes[i];
+		if(pr>=2*nnums)break;
 		int lowoffset=low.divmod(Bigint(pr)).second.lowdigits();
 		int startat;
 		if(lowoffset==0)startat=0;
 		else if((pr-lowoffset)%2==0)startat=(pr-lowoffset)/2;
-		else startat=pr-lowoffset;
+		else startat=(pr-lowoffset+pr)/2;
+		// cerr<<"pr="<<pr<<" lowoffset="<<lowoffset<<" startat="<<startat<<endl;
 		for(int i=startat;i<nnums;i+=pr){ //skips ahead `2*pr` each time (so `pr` array elements)
+			// cerr<<"sieving composite["<<i<<"]: "<<low+2*i<<endl;
+			nleft-=!composite[i];
 			composite[i]=true;
 		}
 	}
+	vector<int> maybeprimes;
+	maybeprimes.reserve(nleft);
+	// cerr<<"Left ("<<nleft<<"):";
+	for(int i=0;i<nnums;i++){
+		if(!composite[i]){
+			// cerr<<' '<<low+2*i;
+			maybeprimes.push_back(i);
+		}
+	}
+	// cerr<<endl;
+
+	while(maybeprimes.size()){
+		int idx=arc4random_uniform(maybeprimes.size());
+		int i=maybeprimes[idx];
+		Bigint bi(low+2*i);
+		if(bailliePSW(bi))return bi;
+		maybeprimes.erase(maybeprimes.begin()+idx);
+	}
+	throw range_error("No primes");
 }
 
-bool bailliePSW(const Bigint&);
+bool strongPseudoPrime2(const Bigint &n){
+	//https://en.wikipedia.org/wiki/Strong_pseudoprime#Formal_definition
+	if(n<2||n.even())return false;
+	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;
+	}
+	if(expmod(Bigint::two,d,n)==1)return true;
+	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;
+	int D;
+	if     ((D=  5),jacobiSymbol(Bigint(D),n)==-1);
+	else if((D= -7),jacobiSymbol(Bigint(D),n)==-1);
+	else if((D=  9),jacobiSymbol(Bigint(D),n)==-1);
+	else if((D=-11),jacobiSymbol(Bigint(D),n)==-1);
+	else {
+		Bigint root(isqrt(n));
+		if(root*root==n)return false; //perfect square
+		for(int i=6;;i++){
+			D=2*i+1;
+			if(i%2==1)D=-D;
+			if(jacobiSymbol(Bigint(D),n)==-1)break;
+		}
+	}
+	//now we have a D
+	int P=1,Q=(1-D)/4;
+	return strongLucasPrime(n,D,P,Q);
+}