aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--bigint.cpp3
-rw-r--r--bigint.h3
-rw-r--r--main.cpp14
-rw-r--r--numalgo.cpp20
-rw-r--r--numalgo.h2
-rw-r--r--primes.cpp70
-rw-r--r--primes.h10
7 files changed, 115 insertions, 7 deletions
diff --git a/bigint.cpp b/bigint.cpp
index 644ffb3..21c4ab2 100644
--- a/bigint.cpp
+++ b/bigint.cpp
@@ -7,9 +7,10 @@
using namespace std;
+Bigint Bigint::mone(-1);
Bigint Bigint::zero(0);
Bigint Bigint::one(1);
-Bigint Bigint::mone(-1);
+Bigint Bigint::two(2);
Bigint::Bigint()
:sign(1){}
diff --git a/bigint.h b/bigint.h
index 213bf5c..76f5415 100644
--- a/bigint.h
+++ b/bigint.h
@@ -89,9 +89,10 @@ public:
digit_t _digit(int idx) const;
+ static Bigint mone;
static Bigint zero;
static Bigint one;
- static Bigint mone;
+ static Bigint two;
};
std::istream& operator>>(std::istream&,Bigint&);
diff --git a/main.cpp b/main.cpp
index 2d652b4..db846a9 100644
--- a/main.cpp
+++ b/main.cpp
@@ -144,10 +144,22 @@ 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;
+}
+
int main(int argc,char **argv){
+ (void)argc;
+ (void)argv;
// biginttest();
// repl(argc,argv);
// performrsa();
// testisqrt(argc,argv);
- fillsmallprimes();
+ // randprime(Bigint(20),Bigint(42));
+ // fermatpseudo();
}
diff --git a/numalgo.cpp b/numalgo.cpp
index f2ebac5..9b1681c 100644
--- a/numalgo.cpp
+++ b/numalgo.cpp
@@ -43,6 +43,26 @@ Bigint expmod(const Bigint &b,const Bigint &e,const Bigint &m){
return res;
}
+int jacobiSymbol(Bigint a,Bigint n){
+ //https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol
+ assert(n>0&&n.odd());
+ int res=1;
+ while(true){
+ 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;
+ while(a.even()){
+ res*=two_over_n;
+ a>>=1;
+ }
+ if(a==1)return res;
+ if(gcd(a,n)!=1)return 0;
+ if(a.divmod(Bigint(4)).second.lowdigits()==3&&nmod8%4==3)res*=-1;
+ swap(a,n);
+ }
+}
+
Bigint isqrt(const Bigint &n){
assert(n>=0);
if(n<=1)return n;
diff --git a/numalgo.h b/numalgo.h
index 997a67c..2304e94 100644
--- a/numalgo.h
+++ b/numalgo.h
@@ -8,6 +8,8 @@ Bigint egcd(const Bigint &a,const Bigint &b,Bigint &x,Bigint &y);
Bigint expmod(const Bigint &base,const Bigint &exponent,const Bigint &modulus);
+int jacobiSymbol(Bigint a,Bigint n);
+
// Returns sqrt(n), rounded down if necessary
Bigint isqrt(const Bigint &n);
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);
+}
diff --git a/primes.h b/primes.h
index 61e5fc7..744fa3e 100644
--- a/primes.h
+++ b/primes.h
@@ -11,9 +11,15 @@ void fillsmallprimes();
//for use in RSA (pass target number of bits of N)
std::pair<Bigint,Bigint> genprimepair(int nbits);
-//finds random in range [low,high]; throws domain_error if no prime found
+//finds random in range [low,high]; throws range_error("No primes") if no prime found
//Will call fillsmallprimes() if not yet done
Bigint randprime(const Bigint &low,const Bigint &high);
+//checks Fermat [pseudo- or actual] primality
+bool strongPseudoPrime2(const Bigint &n);
+
+//checks Lucas [pseudo- or actual] primality
+bool strongLucasPrime(const Bigint &n,int D,int P,int Q);
+
//checks primality
-bool bailliePSW(const Bigint&);
+bool bailliePSW(const Bigint &n);