mpqsPolynom.cc

Go to the documentation of this file.

#include "qsieve-fwd.H"
#include "mpqsPolynom.H"
#include "modulo.H"
#include "mpz_sqrtmod.cc"
#include <cmath>

bool collecting_phase_finished = false;


void CmpqsPolynom::compute_first_polynomial(const mpz_t fuer_kN, const int M)
{
  mpz_set(kN,fuer_kN);
  mpz_div_ui(kN_div2,kN,2); // kN/2 to compute negative modulo values
  
  mpz_mul_ui(D,kN,2); mpz_sqrt(D,D); mpz_div_ui(D,D,2); mpz_div_ui(D,D,M); mpz_sqrt(D,D);
  // D is the initial value for the "optimal" polynomial
  
  // special case:
  //  Problems may occur, if D is (due to small numbers to factorize) a member of the static static factorbase.
  //  Same holds, if D is a dynamic Factor.
  if (mpz_cmp_ui(D,SingleLargePrime_Threshold)<=0)
   {
     if (Factor_Threshold<=1.0)
      {
        collecting_phase_finished = true; // don't collect any dynamic factors!
        cout << "Dynamic Factors are deactivated." << endl;
      }
     else
      cout << "MPQS-polynomial could possibly collide with dynamic factorbase!" << endl;
   }
  if (mpz_cmp_ui(D,StaticFactorbaseSettings::BiggestPrime())<=0)
   {
     cout << "Polynomial collides with static factorbase... enable workaround..." << endl;
     mpz_set_ui(D,StaticFactorbaseSettings::BiggestPrime()+1);
     collecting_phase_finished = true; // don't collect any dynamic factors!
   }

  // now the polynomial computation itself:
  // Let A ~ D^2, D (probab. prime), (D/kN)=1, D=3 (mod 4)
  // for D=3 (mod 4) computation of roots is simpler,
  // but D=1 (mod 2) is also okay.
  // My implementation allows D=1 (mod 2).
  // For D=5 (mod 8) computation of roots is relatively easy, too.
  // and for D=1 (mod 4 resp. 8) one can uses Lucas sequences...

  if (mpz_even_p(D)) mpz_add_ui(D,D,1);

  compute_next_polynomial(); // here the computation will be continued
}

void CmpqsPolynom::compute_next_polynomial(const int step /* =0 */)
{
  mpz_add_ui(D,D,step*4); // "step" helps to distribute polynomial intervals to the clients

  // Let A ~ D^2, D (probab. prime), (D/kN)=1, D=3 (mod 4)
  // for D=3 (mod 4) computation of roots is simpler,
  // but D=1 (mod 2) is also okay.
  // My implementation allows D=1 (mod 2).
  // For D=5 (mod 8) computation of roots is relatively easy, too.
  // and for D=1 (mod 4 resp. 8) one can uses Lucas sequences...

#if 0
  // old method to determine next D
  do mpz_add_ui(D,D,2);
  while ( (!mpz_probab_prime_p(D,probab_prime_checks)) || (mpz_jacobi(D,kN)!=1) );
#else
  // new method to determine next D
  {
    unsigned long int h0=mpz_remainder_ui(D,3*5*7*11*13*17*19*23);
    do
     {
       register unsigned long int h=h0;
       do h+=2; while (!numtheory::coprime(h,3*5*7*11*13*17*19*23));
       mpz_add_ui(D,D,h-h0); h0=h;
       // now D hasn't any small factors
     } while ( (!mpz_probab_prime_p(D,probab_prime_checks)) || (mpz_jacobi(D,kN)!=1) );
    // D is now (probable) prime and no quadratic residue modulo kN
  }
#endif

  mpz_mul(A,D,D); mpz_mul_ui(A2,A,2);

  if (mpz_remainder_ui(D,4)==1) // D=1 (mod 4)
   {
     // h1 = sqrt(kN) (mod D)
     mpz_sqrtmod(h1,kN,D);

     // h0 = h1/kN (mod D) = 1/h1 (mod D)
     mpz_invert(h0,h1,D); mpz_mod(h0,h0,D);
   }
  else // D=3 (mod 4)
   {
     // remark: this case could be handled by the above case, too!
     // But computation is even easier for D=3 (mod 4), because
     // we don't need to invert.

     // h0 = (kN)^((D-3)/4) mod D
     mpz_sub_ui(h0,D,3); mpz_div_ui(h0,h0,4); mpz_powm(h0,kN,h0,D);
  
     // h1 = kN*h0 mod D  (=sqrt(kN) (mod D)) 
     mpz_mul(h1,kN,h0); mpz_mod(h1,h1,D);
   }

  // h2 = (2*h1)^-1 * ((kN-h1^2) /D) mod D
  mpz_mul(h2,h1,h1); mpz_sub(h2,kN,h2);
  
  mpz_t x; mpz_init(x);

  if (!mpz_divisible_p(h2,D))
   {
     MARK;
     mpz_mod(x,h2,D);
     cerr << "Error: Division, unexpected Remainder " << x << " !" << endl;
     exit(1);
   }
      
  mpz_divexact(h2,h2,D);
  mpz_mul_ui(x,h1,2); mpz_invert(x,x,D); mpz_mul(h2,x,h2);
  mpz_mod(h2,h2,D);
  
  // B = h1 + h2*D mod A
  mpz_mul(B,h2,D); mpz_add(B,B,h1); mpz_mod(B,B,A);
  if (mpz_even_p(B)) mpz_sub(B,A,B); 
  
  // C = (B^2-kN) / (4*A)
  mpz_mul(C,B,B); mpz_sub(C,C,kN); mpz_div(C,C,A); mpz_div_ui(C,C,4);

  mpz_clear(x);

#if 0
  // complete, (but unneeded output)
  cout << "Parameters of A*x^2+B*x+C are: " << endl;
  cout << "A=" << A << endl;
  cout << "B=" << B << endl;
  cout << "C=" << C << endl;
  cout << "D=" << D << " with A=D^2" << endl;
#else

#ifdef VERBOSE
  // short output of the essential infos:
  cout << "New polynomial with D=" << D << " computed." << endl;
#endif

#if 0 || defined(PROFILE)
  static unsigned int count = 100;
  if (--count==0) exit(0); // for Profiling-Test
#endif

#endif
  
  mpz_mul_ui(D2_inv_mod_kN,D,2); mpz_invert(D2_inv_mod_kN,D2_inv_mod_kN,kN);
  
  // compute A/D mod kN (which is const. for the selected polynomial)
  mpz_mul_ui(A_div_D_mod_kN,A,2);
  mpz_mul(A_div_D_mod_kN,A_div_D_mod_kN,D2_inv_mod_kN);
  mpz_mod(A_div_D_mod_kN,A_div_D_mod_kN,kN);

  // compute B/2D mod kN (which is const. for the selected polynomial)
  mpz_mul(B_div_2D_mod_kN,B,D2_inv_mod_kN);
  mpz_mod(B_div_2D_mod_kN,B_div_2D_mod_kN,kN);
  
  //testen();
}

void CmpqsPolynom::SanityCheck(const signed int SievePos /* =0 */)
{
  mpz_t x,y;
  mpz_init(x); mpz_init(y);

  mpz_set_si(h2,SievePos);
  
  // Variante 1: ( (2*A*wert+B) / (2*D) )^2 mod kN
  mpz_mul_ui(x,A,2); mpz_mul(x,x,h2); mpz_add(x,x,B);
  //mpz_div(x,x,D); mpz_div_ui(x,x,2);
  mpz_mul(x,x,D2_inv_mod_kN); mpz_mod(x,x,kN);
  mpz_mul(x,x,x); mpz_mod(x,x,kN);
  if (mpz_cmp(x,kN_div2)>0) { mpz_sub(x,kN,x); mpz_neg(x,x); } // compute negative value!
  
  // Variante 2: A*wert^2 + B*wert + C mod kN
  
  mpz_mul(h0,A,h2); mpz_mul(h0,h0,h2);
  mpz_mul(h1,B,h2);
  mpz_add(y,h0,h1); mpz_add(y,y,C); mpz_mod(y,y,kN);
  if (mpz_cmp(y,kN_div2)>0) { mpz_sub(y,kN,y); mpz_neg(y,y); } // compute negative value!
  
  cout << "Test: Q(" << SievePos << ")" << endl;
  cout << " =" << x << endl;
  cout << " =" << y << endl;

  if (mpz_cmp(x,y)!=0)
    {
      cerr << "Error: Values differ!" << endl;
      exit(1);
    }
  mpz_clear(x); mpz_clear(y);
}

void CmpqsPolynom::get_values(const signed int SievePos, mpz_t radix, mpz_t Q) const
{  // Variante 1: ( (2*A*wert+B) / (2*D) )^2 mod kN
  mpz_mul_si(radix,A_div_D_mod_kN,SievePos);
  mpz_add(radix,radix,B_div_2D_mod_kN);
  mpz_mod(radix,radix,kN);
  mpz_mul(Q,radix,radix); mpz_mod(Q,Q,kN);
  if (mpz_cmp(Q,kN_div2)>0) mpz_sub(Q,Q,kN); // compute negative values!
}

double CmpqsPolynom::get_logval(const signed int SievePos) const
{
  mpz_t h;
  mpz_init_set_si(h,SievePos);
  mpz_mul(h,h,A_div_D_mod_kN);
  mpz_add(h,h,B_div_2D_mod_kN);
  mpz_mul(h,h,h);
  mpz_mod(h,h,kN);

  if (mpz_cmp(h,kN_div2)>0) mpz_sub(h,h,kN); // compute negative values!
  const double d=std::log(std::fabs(mpz_get_d(h)));
  mpz_clear(h);
  return d;
}

void CmpqsPolynom::save(ostream &ostr) // for Recovery etc.
{
  // save all data which are important for computing the current polynomial interval
  // D is the value that determines the current polynomial
  ostr << D << endl;
}

void CmpqsPolynom::load(istream &in) // for Recovery etc.
{
  // load all necessary data to compute the current polynomial
  in >> D; in.ignore(1,'\n');
  mpz_sub_ui(D,D,4); // decrease by four
  compute_next_polynomial(); // and compute next polynomial interval
}

void CmpqsPolynom::load_if_available(istream &in) // for Recovery etc.
{
  while (isspace(in.peek())) in.get(); // ignore all whitespaces
  if (in.peek()==EOF)
   {
     cerr << "no mpqs polynomial coefficient found. using start value" << endl;
     compute_next_polynomial(); // using initial values and compute next polynomial interval
   }
  else load(in);
}

---