diff options
| -rw-r--r-- | bigint.cpp | 11 | ||||
| -rw-r--r-- | main.cpp | 12 | ||||
| -rw-r--r-- | numalgo.cpp | 7 | ||||
| -rw-r--r-- | primes.cpp | 74 | ||||
| -rw-r--r-- | primes.h | 2 |
5 files changed, 94 insertions, 12 deletions
@@ -133,8 +133,9 @@ Bigint& Bigint::operator=(slongdigit_t v){ } Bigint& Bigint::operator+=(const Bigint &o){ - if(&o==this){ - return *this=Bigint(*this)+=o; + if(&o==this){ //*this + *this + operator<<=(1); + return *this; } if(sign==1){ if(o.sign==1)add(*this,o); @@ -148,8 +149,10 @@ Bigint& Bigint::operator+=(const Bigint &o){ } Bigint& Bigint::operator-=(const Bigint &o){ - if(&o==this){ - return *this=Bigint(*this)-=o; + if(&o==this){ // *this - *this + sign=1; + digits.clear(); + return *this; } if(sign==1){ if(o.sign==1)subtract(*this,o); @@ -153,6 +153,16 @@ void fermatpseudo(){ cout<<endl; } +void lucaspseudo(){ + fillsmallprimes(); + for(int i=2;i<65000;i++){ + // cerr<<"TRYING "<<i<<endl; + if(!strongLucasPrime(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; @@ -162,4 +172,6 @@ int main(int argc,char **argv){ // testisqrt(argc,argv); // randprime(Bigint(20),Bigint(42)); // fermatpseudo(); + // strongLucasPrime(Bigint(5)); + lucaspseudo(); } diff --git a/numalgo.cpp b/numalgo.cpp index 9b1681c..c3bfa1b 100644 --- a/numalgo.cpp +++ b/numalgo.cpp @@ -8,7 +8,7 @@ Bigint gcd(Bigint a,Bigint b){ while(true){ if(a==0)return b; if(b==0)return a; - if(a>=b)a=a.divmod(b).second; + if(a.compareAbs(b)==1)a=a.divmod(b).second; else b=b.divmod(a).second; } } @@ -45,13 +45,16 @@ Bigint expmod(const Bigint &b,const Bigint &e,const Bigint &m){ int jacobiSymbol(Bigint a,Bigint n){ //https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol + // cerr<<"jacobiSymbol("<<a<<','<<n<<')'<<endl; assert(n>0&&n.odd()); int res=1; while(true){ + // cerr<<" a="<<a<<" n="<<n<<endl; 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; + if(a==0)return 0; while(a.even()){ res*=two_over_n; a>>=1; @@ -84,7 +87,7 @@ Bigint isqrt(const Bigint &n){ if(xnext==x)break; x=xnext; } - cerr<<iter<<endl; + // cerr<<iter<<endl; assert(iter<maxiter); switch((x*x).compare(n)){ case -1:{ @@ -125,9 +125,10 @@ bool strongPseudoPrime2(const Bigint &n){ 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; +bool strongLucasPrime(const Bigint &n){ + //https://en.wikipedia.org/wiki/Lucas_pseudoprime#Implementing_a_Lucas_probable_prime_test + //TODO: use d and s already calculated in strongPseudoPrime2 + if(n.even())return false; int D; if ((D= 5),jacobiSymbol(Bigint(D),n)==-1); else if((D= -7),jacobiSymbol(Bigint(D),n)==-1); @@ -143,6 +144,69 @@ bool bailliePSW(const Bigint &n){ } } //now we have a D - int P=1,Q=(1-D)/4; - return strongLucasPrime(n,D,P,Q); + // cerr<<"D="<<D<<endl; + int P=1,iQ=(1-D)/4; + if(gcd(n,Bigint(P))!=1||gcd(n,Bigint(iQ))!=1)return false; //would be too easy to ignore + + //now begin the actual sequence algorithm +#if 1 + int s=0; + Bigint d(n); + d-=1; + while(d.even()){ + d>>=1; + s++; + } +#else + Bigint d(n); + d+=1; +#endif + // cerr<<"d="<<d<<endl; + vector<bool> dbits=d.bits(); + assert(dbits.size()>0); + assert(dbits[dbits.size()-1]==true); + Bigint U(1),V(P); + Bigint Qk(iQ); + for(int i=dbits.size()-2;i>=0;i--){ + U*=V; + V*=V; + V-=Qk<<1; + Qk*=Qk; + if(dbits[i]){ + Bigint unext(U*P); + unext+=V; + if(unext.odd())unext+=n; + assert(unext.even()); + unext>>=1; + Bigint vnext(U*D); + vnext+=V*P; + if(vnext.odd())vnext+=n; + assert(vnext.even()); + vnext>>=1; + U=unext.divmod(n).second; + V=vnext.divmod(n).second; + Qk=(Qk*iQ).divmod(n).second; + } else { + U=U.divmod(n).second; + V=V.divmod(n).second; + Qk=Qk.divmod(n).second; + } + } + if(U==0)return true; + if(V==0)return true; //r=0 check +#if 1 + for(int r=1;r<s;r++){ + V*=V; + V-=Qk<<1; + V=V.divmod(n).second; + Qk=(Qk*Qk).divmod(n).second; + if(V==0)return true; + } +#endif + return false; +} + +bool bailliePSW(const Bigint &n){ + //https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test#The_test + return strongPseudoPrime2(n)&&strongLucasPrime(n); } @@ -19,7 +19,7 @@ Bigint randprime(const Bigint &low,const Bigint &high); bool strongPseudoPrime2(const Bigint &n); //checks Lucas [pseudo- or actual] primality -bool strongLucasPrime(const Bigint &n,int D,int P,int Q); +bool strongLucasPrime(const Bigint &n); //checks primality bool bailliePSW(const Bigint &n); |