aboutsummaryrefslogtreecommitdiff
path: root/primes.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'primes.cpp')
-rw-r--r--primes.cpp82
1 files changed, 82 insertions, 0 deletions
diff --git a/primes.cpp b/primes.cpp
new file mode 100644
index 0000000..d144e5e
--- /dev/null
+++ b/primes.cpp
@@ -0,0 +1,82 @@
+#include <cstring>
+#include <cmath>
+#include <cassert>
+#include "numalgo.h"
+#include "primes.h"
+
+using namespace std;
+
+vector<int> smallprimes;
+bool smallprimes_inited=false;
+
+void fillsmallprimes(){
+ smallprimes_inited=true;
+ //TODO: reserve expected amount of space in smallprimes
+ smallprimes.push_back(2);
+ const int highbound=65000;
+ bool composite[highbound/2]; //entries for 3, 5, 7, 9, etc.
+ memset(composite,0,highbound/2*sizeof(bool));
+ int roothighbound=sqrt(highbound);
+ for(int i=3;i<=highbound;i+=2){
+ if(composite[i/2-1])continue;
+ smallprimes.push_back(i);
+ if(i>roothighbound)continue;
+ for(int j=i*i;j<=highbound;j+=2*i){
+ composite[j/2-1]=true;
+ }
+ }
+ //for(int p : smallprimes)cerr<<p<<' ';
+ //cerr<<endl;
+}
+
+pair<Bigint,Bigint> genprimepair(int nbits){
+ // for x = nbits/2:
+ // (2^x)^2 = 2^(2x)
+ // (2^x + 2^(x-2))^2 = 2^(2x) + 2^(2x-1) + 2^(2x-4)
+ // ergo: (2^x + lambda*2^(x-2))^2 \in [2^(2x), 2^(2x+1)), for lambda \in [0,1]
+ // To make sure the primes "differ in length by a few digits" [RSA78], we use x1=x-2 in the first
+ // prime and x2-x+2 in the second
+ int x1=nbits/2-2,x2=(nbits+1)/2+2;
+ assert(x1+x2==nbits);
+ return make_pair(
+ randprime(Bigint::one<<x1,(Bigint::one<<x1)+(Bigint::one<<(x1-2))),
+ randprime(Bigint::one<<x2,(Bigint::one<<x2)+(Bigint::one<<(x2-2))));
+}
+
+Bigint randprime(const Bigint &biglow,const Bigint &bighigh){
+ if(!smallprimes_inited)fillsmallprimes();
+
+ assert(bighigh>biglow);
+ static const int maxrangesize=100001;
+ Bigint diff(bighigh-biglow);
+ Bigint low,high; //inclusive
+ if(diff<=maxrangesize){
+ low=biglow;
+ high=bighigh;
+ } else {
+ high=low=cryptrandom_big(diff-maxrangesize);
+ high+=maxrangesize;
+ }
+ if(low.even())low+=1;
+ if(high.even())high-=1;
+ Bigint sizeb(high-low+1);
+ assert(sizeb>=0&&sizeb<=maxrangesize);
+ int nnums=(sizeb.lowdigits()+1)/2;
+
+ 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];
+ 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;
+ for(int i=startat;i<nnums;i+=pr){ //skips ahead `2*pr` each time (so `pr` array elements)
+ composite[i]=true;
+ }
+ }
+}
+
+bool bailliePSW(const Bigint&);