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);
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
00022
00023
00024
00025
00026 if (mpz_cmp_ui(D,SingleLargePrime_Threshold)<=0)
00027 {
00028 if (Factor_Threshold<=1.0)
00029 {
00030 collecting_phase_finished = true;
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;
00041 }
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 if (mpz_even_p(D)) mpz_add_ui(D,D,1);
00052
00053 compute_next_polynomial();
00054 }
00055
00056 void CmpqsPolynom::compute_next_polynomial(const int step )
00057 {
00058 mpz_add_ui(D,D,step*4);
00059
00060
00061
00062
00063
00064
00065
00066
00067 #if 0
00068
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
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
00081 } while ( (!mpz_probab_prime_p(D,probab_prime_checks)) || (mpz_jacobi(D,kN)!=1) );
00082
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)
00089 {
00090
00091 mpz_sqrtmod(h1,kN,D);
00092
00093
00094 mpz_invert(h0,h1,D); mpz_mod(h0,h0,D);
00095 }
00096 else
00097 {
00098
00099
00100
00101
00102
00103 mpz_sub_ui(h0,D,3); mpz_div_ui(h0,h0,4); mpz_powm(h0,kN,h0,D);
00104
00105
00106 mpz_mul(h1,kN,h0); mpz_mod(h1,h1,D);
00107 }
00108
00109
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
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
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
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
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);
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
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
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
00168 }
00169
00170 void CmpqsPolynom::SanityCheck(const signed int SievePos )
00171 {
00172 mpz_t x,y;
00173 mpz_init(x); mpz_init(y);
00174
00175 mpz_set_si(h2,SievePos);
00176
00177
00178 mpz_mul_ui(x,A,2); mpz_mul(x,x,h2); mpz_add(x,x,B);
00179
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); }
00183
00184
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); }
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 {
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);
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);
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)
00228 {
00229
00230
00231 ostr << D << endl;
00232 }
00233
00234 void CmpqsPolynom::load(istream &in)
00235 {
00236
00237 in >> D; in.ignore(1,'\n');
00238 mpz_sub_ui(D,D,4);
00239 compute_next_polynomial();
00240 }
00241
00242 void CmpqsPolynom::load_if_available(istream &in)
00243 {
00244 while (isspace(in.peek())) in.get();
00245 if (in.peek()==EOF)
00246 {
00247 cerr << "no mpqs polynomial coefficient found. using start value" << endl;
00248 compute_next_polynomial();
00249 }
00250 else load(in);
00251 }