diff options
| -rw-r--r-- | bigint.cpp | 3 | ||||
| -rw-r--r-- | bigint.h | 3 | ||||
| -rw-r--r-- | main.cpp | 14 | ||||
| -rw-r--r-- | numalgo.cpp | 20 | ||||
| -rw-r--r-- | numalgo.h | 2 | ||||
| -rw-r--r-- | primes.cpp | 70 | ||||
| -rw-r--r-- | primes.h | 10 |
7 files changed, 115 insertions, 7 deletions
@@ -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){} @@ -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&); @@ -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; @@ -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); @@ -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); +} @@ -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); |