aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--bigint.cpp11
-rw-r--r--main.cpp12
-rw-r--r--numalgo.cpp7
-rw-r--r--primes.cpp74
-rw-r--r--primes.h2
5 files changed, 94 insertions, 12 deletions
diff --git a/bigint.cpp b/bigint.cpp
index 21c4ab2..9fdafe0 100644
--- a/bigint.cpp
+++ b/bigint.cpp
@@ -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);
diff --git a/main.cpp b/main.cpp
index db846a9..43d2de7 100644
--- a/main.cpp
+++ b/main.cpp
@@ -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:{
diff --git a/primes.cpp b/primes.cpp
index 135e7bb..c9d25c0 100644
--- a/primes.cpp
+++ b/primes.cpp
@@ -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);
}
diff --git a/primes.h b/primes.h
index 744fa3e..eb817b1 100644
--- a/primes.h
+++ b/primes.h
@@ -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);