#include "polphi_template.H"
#include "my_mpz_powm_ui.cc"
/*! @file
* @brief
* contains implementation classes for instantiating polphi-template to pollard phi method
*/
class CRingPhi
{
private:
mpz_t a;
public:
static const std::string name;
static const int SieveIntervalSize=2*3*5*7*11*13*17;
CRingPhi() { mpz_init(a); }
CRingPhi(const CRingPhi &x) { mpz_init_set(a,x.a); }
~CRingPhi() { mpz_clear(a); }
bool set_startvalue(const unsigned int seed=1) { mpz_set_ui(a,seed+1); return true; }
void set(const CRingPhi &x) { mpz_set(a,x.a); }
void pow_mod(const unsigned int i, const mpz_t M) { mpz_powm_ui(a,a,i,M); }
void fast_pow_mod(const unsigned int i, const mpz_t M) { my_mpz_powm_ui(a,i,M); }
void test_gcd(mpz_t x, const mpz_t M) { mpz_sub_ui(x,a,1); mpz_gcd(x,x,M); }
friend class CRingPhiPhase2;
};
const std::string CRingPhi::name = "phi";
class CRingPhiPhase2 : private ForbidAssignment
{
private:
CRingPhi a;
mpz_t b_inv, b, d, d_inv, b_plus_b_inv, d_plus_d_inv, prev_b_plus_b_inv,
temp1, temp2;
public:
explicit CRingPhiPhase2(const CRingPhi &_a) : a(_a)
{
// "a" contains the starting value (base) for Phase 2
/* rough description:
utilize the identity:
(a^(i+d) -1) * (a^(i-d) -1) = a^i * ( (a^i +a^(-i)) - (a^d +a^(-d)) )
1. The multiplier "a^i" can be ignored, as it has no influence on the gcd
2. D will be increased in steps of i.
-> precalculate (a^d +a^(-d)) for 0<d<D,
-> furthermore: per interval determine/update (a^i + a^(-i))
3. a^(i+D) can be computed by using a^(i+D)=a^i * a^D;
a^(-(i+D)) can be computed analogously, because a^(-i)=(a^(-1))^i.
-> two multiplications per interval to determine a^(i+D) + a^(-(i+D)).
Bemerkung I: (remark I; not very important)
Die hierauf aufbauende Paarbildungsmethode konstruiert nicht
die durch andere Methoden erreichbaren Primzahlpaare, sie ist
aber dafür vergleichsweise einfach zu implementieren und sie
kommt mit weniger Speicherplatz als die improved-standard-continuation
aus. Sobald mehrere Primzahlpaare innerhalb eines Intervalls
auftauchen, ist sie schneller als die improved standard continuation.
remark II: (ad 3.)
a^i+a^(-i) can be computed efficiently via Lucas sequences
because of x^2-x(a^1+a^(-1))+1=(x-a)*(x-a^(-1)):
a^i+a^(-i)=Vi(a^1+a^(-1))=:V[i];
a^(i+D)+a^(-(i+D)))=V[i+D]=V[i]*V[D]-V[i-D].
The number of required multiplications per interval-step can be reduced from 2 to 1.
This is another way leading (conceptually) to the same continuation as it is
described in the paper of P.L. Montgomery.
(cf. Peter L. Montgomery, "Speeding the Pollard and Elliptic Curve
Methods of Factorization", MathComp Vol 48 (177), Jan 1987, p.243-264;
especially section 4, 5 and 9).
*/
mpz_init(b); mpz_init(b_inv); mpz_init(d); mpz_init(d_inv); mpz_init(b_plus_b_inv);
mpz_init(d_plus_d_inv); mpz_init(prev_b_plus_b_inv);
mpz_init(temp1); mpz_init(temp2);
// b,b_inv, d,d_inv don't get updated anymore in the following,
// because the associated sums will be computed via Lucas sequences.
// Precomputing most common used values:
mpz_set(b,a.a); mpz_invert(b_inv,b,n); mpz_mod(b_inv,b_inv,n);
mpz_powm_ui(d,b,2,n); mpz_powm_ui(d_inv,b_inv,2,n);
mpz_add(d_plus_d_inv,d,d_inv); mpz_mod(d_plus_d_inv,d_plus_d_inv,n);
mpz_add(b_plus_b_inv,b,b_inv); mpz_mod(b_plus_b_inv,b_plus_b_inv,n);
mpz_set(prev_b_plus_b_inv,b_plus_b_inv);
}
~CRingPhiPhase2()
{
mpz_clear(b); mpz_clear(b_inv); mpz_clear(d); mpz_clear(d_inv); mpz_clear(b_plus_b_inv);
mpz_clear(d_plus_d_inv); mpz_clear(prev_b_plus_b_inv);
mpz_clear(temp1); mpz_clear(temp2);
}
void get_polynomdef_point(mpz_t x)
{
mpz_set(x,b_plus_b_inv);
}
void calc_polynomdef_next_point()
{
mpz_set(temp1,prev_b_plus_b_inv); mpz_set(prev_b_plus_b_inv,b_plus_b_inv);
mpz_mul(b_plus_b_inv,b_plus_b_inv,d_plus_d_inv);
mpz_sub(b_plus_b_inv,b_plus_b_inv,temp1);
mpz_mod(b_plus_b_inv,b_plus_b_inv,n);
}
void calc_EvalStartingPoint(const int D, const double ii)
{
mpz_powm_ui(d,a.a,D,n); mpz_invert(d_inv,d,n); mpz_mod(d_inv,d_inv,n);
mpz_add(d_plus_d_inv,d,d_inv); mpz_mod(d_plus_d_inv,d_plus_d_inv,n);
mpz_set_d(temp1,ii); mpz_powm(b,a.a,temp1,n);
mpz_invert(b_inv,b,n); mpz_mod(b_inv,b_inv,n);
mpz_add(b_plus_b_inv,b,b_inv); mpz_mod(b_plus_b_inv,b_plus_b_inv,n);
mpz_mul(temp1,b,d_inv); mpz_mod(temp1,temp1,n);
mpz_mul(temp2,b_inv,d); mpz_mod(temp2,temp2,n);
mpz_add(prev_b_plus_b_inv,temp1,temp2); mpz_mod(prev_b_plus_b_inv,prev_b_plus_b_inv,n);
}
void get_point_and_calc_next_point(mpz_t x)
{
mpz_set(x,prev_b_plus_b_inv); mpz_set(prev_b_plus_b_inv,b_plus_b_inv);
mpz_mul(b_plus_b_inv,b_plus_b_inv,d_plus_d_inv);
mpz_sub(b_plus_b_inv,b_plus_b_inv,x);
mpz_mod(b_plus_b_inv,b_plus_b_inv,n);
mpz_set(x,prev_b_plus_b_inv);
}
};
// explicit instantiation (not necessary)
template void polphi_template<CRingPhi,CRingPhiPhase2>(const int phase1, const double phase2);
#define pollard_phi(phase1,phase2) polphi_template<CRingPhi,CRingPhiPhase2>(phase1,phase2)