Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
mersenne first factor
#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <cstdlib>
#include <stdint.h>
int64_t Euclidean_Diophantine( int64_t coeff_a, int64_t coeff_b)
{
// return y in Diophantine equation coeff_a x + coeff_b y = 1
int64_t k=1;
std::vector<int64_t> div_t;
std::vector<int64_t> rem_t;
std::vector<int64_t> coeff_t;
div_t.push_back(coeff_a);
rem_t.push_back(coeff_b);
coeff_t.push_back((int64_t)0);
div_t.push_back((int64_t)div_t[0]/rem_t[0]);
rem_t.push_back((int64_t)div_t[0]%rem_t[0]);
coeff_t.push_back((int64_t)0);
while (rem_t[k]>1)
{
k=k+1;
div_t.push_back((int64_t)rem_t[k-2]/rem_t[k-1]);
rem_t.push_back((int64_t)rem_t[k-2]%rem_t[k-1]);
coeff_t.push_back((int64_t)0);
}
k=k-1;
coeff_t[k]=-div_t[k+1];
if (k>0)
coeff_t[k-1]=(int64_t)1;
while (k > 1)
{
k=k-1;
coeff_t[k-1]=coeff_t[k+1];
coeff_t[k]+=(int64_t)(coeff_t[k+1]*(-div_t[k+1]));
}
if (k==1)
return (int64_t)(coeff_t[k-1]+coeff_t[k]*(-div_t[k]));
else
return (int64_t)(coeff_t[0]);
}
uint64_t XY_modulo_D_LN(uint64_t X,uint64_t Y,uint64_t D,uint64_t D_max)
{
//D>D_max tested for D<2^54
uint64_t mod_D;
if (X<=D_max && Y<=D_max )
{
mod_D=(X*Y)%D;
}
else if (X<=D_max || Y<=D_max )
{
uint64_t N_t;
uint64_t RD=(uint64_t)D%D_max;
uint64_t KD=(uint64_t)D/D_max;
if (Y<=D_max)
{
uint64_t X1=Y;
Y=X;
X=X1;
}
mod_D=(X*(Y%D_max))%D_max;
N_t=X*(Y/D_max)+(X*(Y%D_max))/D_max;
mod_D=(mod_D+(D_max*(N_t%KD))%D)%D;
mod_D=(mod_D+D-((N_t/KD)*RD)%D)%D;
}
else
{
uint64_t RD=(uint64_t)D%D_max;
uint64_t KD=(uint64_t)D/D_max;
uint64_t Rx=(uint64_t)X%D_max;
uint64_t Ry=(uint64_t)Y%D_max;
uint64_t Kx=(uint64_t)X/D_max;
uint64_t Ky=(uint64_t)Y/D_max;
uint64_t mod_D_t,N_t,R_t;
mod_D=(Rx*Ry)%D;
mod_D_t=((D_max*(D_max%KD))%D)%D;
mod_D_t=(mod_D_t+D-((D_max/KD)*RD)%D)%D;
N_t=mod_D_t;
R_t=Kx;
if (N_t>=(D_max/R_t))
{
mod_D_t=(R_t*(N_t%D_max))%D_max;
N_t=R_t*(N_t/D_max)+(R_t*(N_t%D_max))/D_max;
mod_D_t=(mod_D_t+(D_max*(N_t%KD))%D)%D;
mod_D_t=(mod_D_t+D-((N_t/KD)*RD)%D)%D;
}
else
mod_D_t=mod_D_t*R_t;
N_t=mod_D_t;
R_t=Ky;
if (N_t>=(D_max/R_t))
{
mod_D_t=(R_t*(N_t%D_max))%D_max;
N_t=R_t*(N_t/D_max)+(R_t*(N_t%D_max))/D_max;
mod_D_t=(mod_D_t+(D_max*(N_t%KD))%D)%D;
mod_D_t=(mod_D_t+D-((N_t/KD)*RD)%D)%D;
}
else
mod_D_t=mod_D_t*R_t;
N_t=mod_D+mod_D_t;
mod_D=N_t;
if (N_t>=D_max)
{
mod_D=(N_t)%D_max;
mod_D=(mod_D+(D_max*((N_t/D_max)%KD))%D)%D;
mod_D=(mod_D+D-(((N_t/D_max)/KD)*RD)%D)%D;
}
if (Rx>=KD)
{
mod_D_t=((D_max*(Rx%KD))%D)%D;
mod_D_t=(mod_D_t+D-((Rx/KD)*RD)%D)%D;
}
else
mod_D_t=(Rx*D_max)%D;
N_t=mod_D_t;
R_t=Ky;
if (N_t>=(D_max/R_t))
{
mod_D_t=(R_t*(N_t%D_max))%D_max;
N_t=R_t*(N_t/D_max)+(R_t*(N_t%D_max))/D_max;
mod_D_t=(mod_D_t+(D_max*(N_t%KD))%D)%D;
mod_D_t=(mod_D_t+D-((N_t/KD)*RD)%D)%D;
}
else
mod_D_t=mod_D_t*R_t;
N_t=mod_D+mod_D_t;
mod_D=N_t;
if (N_t>=D_max)
{
mod_D=(N_t)%D_max;
mod_D=(mod_D+(D_max*((N_t/D_max)%KD))%D)%D;
mod_D=(mod_D+D-(((N_t/D_max)/KD)*RD)%D)%D;
}
if (Ry>=KD)
{
mod_D_t=((D_max*(Ry%KD))%D)%D;
mod_D_t=(mod_D_t+D-((Ry/KD)*RD)%D)%D;
}
else
mod_D_t=(Ry*D_max)%D;
N_t=mod_D_t;
R_t=Kx;
if (N_t>=(D_max/R_t))
{
mod_D_t=(R_t*(N_t%D_max))%D_max;
N_t=R_t*(N_t/D_max)+(R_t*(N_t%D_max))/D_max;
mod_D_t=(mod_D_t+(D_max*(N_t%KD))%D)%D;
mod_D_t=(mod_D_t+D-((N_t/KD)*RD)%D)%D;
}
else
mod_D_t=mod_D_t*R_t;
N_t=mod_D+mod_D_t;
mod_D=N_t;
if (N_t>=D_max)
{
mod_D=(N_t)%D_max;
mod_D=(mod_D+(D_max*((N_t/D_max)%KD))%D)%D;
mod_D=(mod_D+D-(((N_t/D_max)/KD)*RD)%D)%D;
}
}
return mod_D;
}
bool is_2pow_onemod(uint64_t P,uint64_t D,uint64_t D_max)
{
// return true if 2^P = 1 (mod D)
uint64_t x=2;
uint64_t y=1;
if (D<=D_max)
{
while (P > (uint64_t)1)
{
if (P%2 == 1)
y = (x * y) %D;
x = (x * x) %D;
P >>= 1;
}
return (((x*y)%D) == 1);
}
else
{
while (P > (uint64_t)1)
{
if (P%2 == 1)
y=XY_modulo_D_LN(x,y,D,D_max);
x=XY_modulo_D_LN(x,x,D,D_max);
P >>= 1;
}
return (XY_modulo_D_LN(x,y,D,D_max) == (uint64_t)1);
}
}
uint64_t mersenne_first_factor(uint64_t P)
{
// P>61 returns, if it exists, the first factor D of (2^P - 1) less than 1+dim_step*8*P
int64_t bW=(int64_t)(8*P);
// set the maximum value for factor D==1+dim_step*bW tested for D<=2^46
int64_t dim_step=(int64_t)2097152;
// set the maximum value small primes for sieve
int64_t dim_primes=std::min(bW,(int64_t)524288);
// set the maximum value for modulo function XY_modulo_D_LN
uint64_t D_max_mod=(uint64_t)1<<30;
uint64_t D1mod8=(uint64_t)1;
uint64_t D7mod8=(uint64_t)1+2*P;
if ((P/(uint64_t)6)%2==0)
D7mod8+=(uint64_t)4*P;
// vector for find small primes for sieve
std::vector<char> Primes_s(dim_primes + 1,true);
// vector for find factor primes 1 mod 8 p_k= 1+8*P*k for k>0
std::vector<char> Primes_1mod8(dim_step + 1,true);
// vector for find factor 7 mod 8 p_k= 1+2*p+8*P*k or p_k= 1+6*p+8*P*k for k>=0
std::vector<char> Primes_7mod8(dim_step + 1,true);
int64_t r_D,r_t,r_t1,C1,C2,m, mmin;
// find small primes for sieve
for (int64_t p=3;p*p<=dim_primes;p+=2)
if (Primes_s[p])
for (m=p*p;m<=dim_primes;m+=p)
Primes_s[m]=false;
if ((int64_t)P<=dim_primes)
Primes_s[(int64_t)P]=false;
if ((int64_t)D7mod8<=dim_primes)
Primes_s[(int64_t)D7mod8]=false;
r_D=(int64_t)D7mod8;
for (int64_t p=3; p<dim_primes;p+=2)
{
if (Primes_s[p])
{
r_t1=Euclidean_Diophantine(bW,p);
// delete multiples of small primes in Primes_1mod8
r_t=(int64_t)r_t1%bW;
if (r_t>0)
r_t-=bW;
C2=(int64_t)((-bW+p)*r_t/bW);
C1=r_t-bW+p;
mmin=(int64_t)(bW+C1+C2);
for (m =mmin; m<=dim_step && m>=0; m += p)
Primes_1mod8[m] = false;
// delete multiples of small primes in Primes_7mod8
r_t=(int64_t)(r_t1*r_D)%bW;
if (r_t>0)
r_t-=bW;
C2=(int64_t)((-bW+p)*r_t/bW);
C1=r_t-bW+p;
mmin=(int64_t)(bW+C1+C2);
for (m =mmin; m<=dim_step && m>=0; m += p)
Primes_7mod8[m] = false;
}
}
Primes_1mod8[0]=false;
Primes_7mod8[0]=true;
uint64_t D=(uint64_t)0;
for (int64_t k=0; k<=dim_step; D1mod8+=(uint64_t)8*P, D7mod8+=(uint64_t)8*P, k++)
{
if (Primes_7mod8[k])
{
if (is_2pow_onemod(P,D7mod8,D_max_mod))
{
D=D7mod8;
break;
}
}
if (Primes_1mod8[k])
{
if (is_2pow_onemod(P,D1mod8,D_max_mod))
{
D=D1mod8;
break;
}
}
}
return D;
}
int main()
{
std::cout <<"factor: "<< mersenne_first_factor(668267879)<< std::endl;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment