pollard_phi.cc

Go to the documentation of this file.

#include "polphi_template.H"
#include "my_mpz_powm_ui.cc"

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)

---