mpqsPolynom.cc

Go to the documentation of this file.
00001 
00006 #include "qsieve-fwd.H"
00007 #include "mpqsPolynom.H"
00008 #include "modulo.H"
00009 #include "mpz_sqrtmod.cc"
00010 #include <cmath>
00011 
00012 bool collecting_phase_finished = false;
00013 
00014 
00015 void CmpqsPolynom::compute_first_polynomial(const mpz_t fuer_kN, const int M)
00016 {
00017   mpz_set(kN,fuer_kN);
00018   mpz_div_ui(kN_div2,kN,2); // kN/2 to compute negative modulo values
00019   
00020   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);
00021   // D is the initial value for the "optimal" polynomial
00022   
00023   // special case:
00024   //  Problems may occur, if D is (due to small numbers to factorize) a member of the static static factorbase.
00025   //  Same holds, if D is a dynamic Factor.
00026   if (mpz_cmp_ui(D,SingleLargePrime_Threshold)<=0)
00027    {
00028      if (Factor_Threshold<=1.0)
00029       {
00030         collecting_phase_finished = true; // don't collect any dynamic factors!
00031         cout << "Dynamic Factors are deactivated." << endl;
00032       }
00033      else
00034       cout << "MPQS-polynomial could possibly collide with dynamic factorbase!" << endl;
00035    }
00036   if (mpz_cmp_ui(D,StaticFactorbaseSettings::BiggestPrime())<=0)
00037    {
00038      cout << "Polynomial collides with static factorbase... enable workaround..." << endl;
00039      mpz_set_ui(D,StaticFactorbaseSettings::BiggestPrime()+1);
00040      collecting_phase_finished = true; // don't collect any dynamic factors!
00041    }
00042 
00043   // now the polynomial computation itself:
00044   // Let A ~ D^2, D (probab. prime), (D/kN)=1, D=3 (mod 4)
00045   // for D=3 (mod 4) computation of roots is simpler,
00046   // but D=1 (mod 2) is also okay.
00047   // My implementation allows D=1 (mod 2).
00048   // For D=5 (mod 8) computation of roots is relatively easy, too.
00049   // and for D=1 (mod 4 resp. 8) one can uses Lucas sequences...
00050 
00051   if (mpz_even_p(D)) mpz_add_ui(D,D,1);
00052 
00053   compute_next_polynomial(); // here the computation will be continued
00054 }
00055 
00056 void CmpqsPolynom::compute_next_polynomial(const int step /* =0 */)
00057 {
00058   mpz_add_ui(D,D,step*4); // "step" helps to distribute polynomial intervals to the clients
00059 
00060   // Let A ~ D^2, D (probab. prime), (D/kN)=1, D=3 (mod 4)
00061   // for D=3 (mod 4) computation of roots is simpler,
00062   // but D=1 (mod 2) is also okay.
00063   // My implementation allows D=1 (mod 2).
00064   // For D=5 (mod 8) computation of roots is relatively easy, too.
00065   // and for D=1 (mod 4 resp. 8) one can uses Lucas sequences...
00066 
00067 #if 0
00068   // old method to determine next D
00069   do mpz_add_ui(D,D,2);
00070   while ( (!mpz_probab_prime_p(D,probab_prime_checks)) || (mpz_jacobi(D,kN)!=1) );
00071 #else
00072   // new method to determine next D
00073   {
00074     unsigned long int h0=mpz_remainder_ui(D,3*5*7*11*13*17*19*23);
00075     do
00076      {
00077        register unsigned long int h=h0;
00078        do h+=2; while (!numtheory::coprime(h,3*5*7*11*13*17*19*23));
00079        mpz_add_ui(D,D,h-h0); h0=h;
00080        // now D hasn't any small factors
00081      } while ( (!mpz_probab_prime_p(D,probab_prime_checks)) || (mpz_jacobi(D,kN)!=1) );
00082     // D is now (probable) prime and no quadratic residue modulo kN
00083   }
00084 #endif
00085 
00086   mpz_mul(A,D,D); mpz_mul_ui(A2,A,2);
00087 
00088   if (mpz_remainder_ui(D,4)==1) // D=1 (mod 4)
00089    {
00090      // h1 = sqrt(kN) (mod D)
00091      mpz_sqrtmod(h1,kN,D);
00092 
00093      // h0 = h1/kN (mod D) = 1/h1 (mod D)
00094      mpz_invert(h0,h1,D); mpz_mod(h0,h0,D);
00095    }
00096   else // D=3 (mod 4)
00097    {
00098      // remark: this case could be handled by the above case, too!
00099      // But computation is even easier for D=3 (mod 4), because
00100      // we don't need to invert.
00101 
00102      // h0 = (kN)^((D-3)/4) mod D
00103      mpz_sub_ui(h0,D,3); mpz_div_ui(h0,h0,4); mpz_powm(h0,kN,h0,D);
00104   
00105      // h1 = kN*h0 mod D  (=sqrt(kN) (mod D)) 
00106      mpz_mul(h1,kN,h0); mpz_mod(h1,h1,D);
00107    }
00108 
00109   // h2 = (2*h1)^-1 * ((kN-h1^2) /D) mod D
00110   mpz_mul(h2,h1,h1); mpz_sub(h2,kN,h2);
00111   
00112   mpz_t x; mpz_init(x);
00113 
00114   if (!mpz_divisible_p(h2,D))
00115    {
00116      MARK;
00117      mpz_mod(x,h2,D);
00118      cerr << "Error: Division, unexpected Remainder " << x << " !" << endl;
00119      exit(1);
00120    }
00121       
00122   mpz_divexact(h2,h2,D);
00123   mpz_mul_ui(x,h1,2); mpz_invert(x,x,D); mpz_mul(h2,x,h2);
00124   mpz_mod(h2,h2,D);
00125   
00126   // B = h1 + h2*D mod A
00127   mpz_mul(B,h2,D); mpz_add(B,B,h1); mpz_mod(B,B,A);
00128   if (mpz_even_p(B)) mpz_sub(B,A,B); 
00129   
00130   // C = (B^2-kN) / (4*A)
00131   mpz_mul(C,B,B); mpz_sub(C,C,kN); mpz_div(C,C,A); mpz_div_ui(C,C,4);
00132 
00133   mpz_clear(x);
00134 
00135 #if 0
00136   // complete, (but unneeded output)
00137   cout << "Parameters of A*x^2+B*x+C are: " << endl;
00138   cout << "A=" << A << endl;
00139   cout << "B=" << B << endl;
00140   cout << "C=" << C << endl;
00141   cout << "D=" << D << " with A=D^2" << endl;
00142 #else
00143 
00144 #ifdef VERBOSE
00145   // short output of the essential infos:
00146   cout << "New polynomial with D=" << D << " computed." << endl;
00147 #endif
00148 
00149 #if 0 || defined(PROFILE)
00150   static unsigned int count = 100;
00151   if (--count==0) exit(0); // for Profiling-Test
00152 #endif
00153 
00154 #endif
00155   
00156   mpz_mul_ui(D2_inv_mod_kN,D,2); mpz_invert(D2_inv_mod_kN,D2_inv_mod_kN,kN);
00157   
00158   // compute A/D mod kN (which is const. for the selected polynomial)
00159   mpz_mul_ui(A_div_D_mod_kN,A,2);
00160   mpz_mul(A_div_D_mod_kN,A_div_D_mod_kN,D2_inv_mod_kN);
00161   mpz_mod(A_div_D_mod_kN,A_div_D_mod_kN,kN);
00162 
00163   // compute B/2D mod kN (which is const. for the selected polynomial)
00164   mpz_mul(B_div_2D_mod_kN,B,D2_inv_mod_kN);
00165   mpz_mod(B_div_2D_mod_kN,B_div_2D_mod_kN,kN);
00166   
00167   //testen();
00168 }
00169 
00170 void CmpqsPolynom::SanityCheck(const signed int SievePos /* =0 */)
00171 {
00172   mpz_t x,y;
00173   mpz_init(x); mpz_init(y);
00174 
00175   mpz_set_si(h2,SievePos);
00176   
00177   // Variante 1: ( (2*A*wert+B) / (2*D) )^2 mod kN
00178   mpz_mul_ui(x,A,2); mpz_mul(x,x,h2); mpz_add(x,x,B);
00179   //mpz_div(x,x,D); mpz_div_ui(x,x,2);
00180   mpz_mul(x,x,D2_inv_mod_kN); mpz_mod(x,x,kN);
00181   mpz_mul(x,x,x); mpz_mod(x,x,kN);
00182   if (mpz_cmp(x,kN_div2)>0) { mpz_sub(x,kN,x); mpz_neg(x,x); } // compute negative value!
00183   
00184   // Variante 2: A*wert^2 + B*wert + C mod kN
00185   
00186   mpz_mul(h0,A,h2); mpz_mul(h0,h0,h2);
00187   mpz_mul(h1,B,h2);
00188   mpz_add(y,h0,h1); mpz_add(y,y,C); mpz_mod(y,y,kN);
00189   if (mpz_cmp(y,kN_div2)>0) { mpz_sub(y,kN,y); mpz_neg(y,y); } // compute negative value!
00190   
00191   cout << "Test: Q(" << SievePos << ")" << endl;
00192   cout << " =" << x << endl;
00193   cout << " =" << y << endl;
00194 
00195   if (mpz_cmp(x,y)!=0)
00196     {
00197       cerr << "Error: Values differ!" << endl;
00198       exit(1);
00199     }
00200   mpz_clear(x); mpz_clear(y);
00201 }
00202 
00203 void CmpqsPolynom::get_values(const signed int SievePos, mpz_t radix, mpz_t Q) const
00204 {  // Variante 1: ( (2*A*wert+B) / (2*D) )^2 mod kN
00205   mpz_mul_si(radix,A_div_D_mod_kN,SievePos);
00206   mpz_add(radix,radix,B_div_2D_mod_kN);
00207   mpz_mod(radix,radix,kN);
00208   mpz_mul(Q,radix,radix); mpz_mod(Q,Q,kN);
00209   if (mpz_cmp(Q,kN_div2)>0) mpz_sub(Q,Q,kN); // compute negative values!
00210 }
00211 
00212 double CmpqsPolynom::get_logval(const signed int SievePos) const
00213 {
00214   mpz_t h;
00215   mpz_init_set_si(h,SievePos);
00216   mpz_mul(h,h,A_div_D_mod_kN);
00217   mpz_add(h,h,B_div_2D_mod_kN);
00218   mpz_mul(h,h,h);
00219   mpz_mod(h,h,kN);
00220 
00221   if (mpz_cmp(h,kN_div2)>0) mpz_sub(h,h,kN); // compute negative values!
00222   const double d=std::log(std::fabs(mpz_get_d(h)));
00223   mpz_clear(h);
00224   return d;
00225 }
00226 
00227 void CmpqsPolynom::save(ostream &ostr) // for Recovery etc.
00228 {
00229   // save all data which are important for computing the current polynomial interval
00230   // D is the value that determines the current polynomial
00231   ostr << D << endl;
00232 }
00233 
00234 void CmpqsPolynom::load(istream &in) // for Recovery etc.
00235 {
00236   // load all necessary data to compute the current polynomial
00237   in >> D; in.ignore(1,'\n');
00238   mpz_sub_ui(D,D,4); // decrease by four
00239   compute_next_polynomial(); // and compute next polynomial interval
00240 }
00241 
00242 void CmpqsPolynom::load_if_available(istream &in) // for Recovery etc.
00243 {
00244   while (isspace(in.peek())) in.get(); // ignore all whitespaces
00245   if (in.peek()==EOF)
00246    {
00247      cerr << "no mpqs polynomial coefficient found. using start value" << endl;
00248      compute_next_polynomial(); // using initial values and compute next polynomial interval
00249    }
00250   else load(in);
00251 }

Generated on Wed Nov 7 23:29:25 2007 for Qsieve by  doxygen 1.5.4