diff options
Diffstat (limited to 'primes.cpp')
-rw-r--r-- | primes.cpp | 82 |
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&); |