00001 #ifndef polphi_template_HEADER
00002 #define polphi_template_HEADER
00003
00010 #include <vector>
00011 #include <sstream>
00012
00013 extern void get_fft_parameter(unsigned int &pg_returnvalue, unsigned int &pms_returnvalue,
00014 const double phase2, const mpz_t N);
00015
00016 template <class CRing, class CRingPhase2>
00017 void polphi_template(const int phase1, const double phase2)
00018 {
00019
00020
00021
00022
00023
00024
00025
00026
00027 using numtheory::is_prime;
00028 using numtheory::gcd;
00029
00030 mpz_t x,y;
00031 mpz_init(x); mpz_init(y);
00032
00033 CRing a,b;
00034
00035 std::vector<int> KnownIndices;
00036 int PreviousIndex = 0;
00037 unsigned int polphi_seed = 0;
00038 const int i_Start = 2;
00039 int i;
00040 short int d;
00041
00042 restart_phi0_with_new_seed:
00043 KnownIndices.clear(); PreviousIndex=0;
00044 ++polphi_seed;
00045
00046
00047
00048
00049
00050 cout << "--------------------------------------------------" << endl;
00051 cout << CRing::name << "-phase 0" << endl;
00052 if (!a.set_startvalue(polphi_seed)) goto done;
00053 for (i=i_Start; i<10000; i++)
00054 {
00055 a.pow_mod(i,n);
00056 retry:
00057 a.test_gcd(x,n);
00058 if (mpz_cmp_ui(x,1)!=0)
00059 {
00060 cout << "factor found in phase 0: i=" << i << endl;
00061 if (mpz_cmp(x,n)==0)
00062 {
00063 cout << "factor is trivial." << endl;
00064
00065 int h=i; int j=2;
00066 while (h>j)
00067 {
00068 while(h>j && h%j==0) h/=j;
00069 ++j;
00070 }
00071
00072 if (KnownIndices.empty()) KnownIndices.push_back(h);
00073 else
00074 if (PreviousIndex!=i || KnownIndices.back()!=h) KnownIndices.push_back(h);
00075 else
00076 {
00077
00078 if (KnownIndices.back()!=i/h) KnownIndices.pop_back();
00079
00080 KnownIndices.push_back(i/h);
00081 }
00082 PreviousIndex=i;
00083
00084 #ifdef VERBOSE_INFO
00085 cout << "KnownIndices:";
00086 for (std::vector<int>::const_iterator p=KnownIndices.begin(); p!=KnownIndices.end(); ++p) cout << " " << *p;
00087 cout << endl;
00088 #endif
00089
00090 if (KnownIndices.empty() || KnownIndices.back()==1)
00091 {
00092 #ifdef VERBOSE_INFO
00093 MARK; cout << "Giving up for " << n << endl;
00094
00095 #endif
00096 goto restart_phi0_with_new_seed;
00097 }
00098
00099 if (KnownIndices.size()<100)
00100 {
00101 a.set_startvalue(polphi_seed);
00102 for (std::vector<int>::const_iterator p=KnownIndices.begin(); p!=KnownIndices.end(); ++p) a.pow_mod(*p,n);
00103 cout << "Trying again..." << endl;
00104 i=i_Start; goto retry;
00105 }
00106 #ifdef VERBOSE_WARN
00107 MARK; cout << "Aborting for " << n << endl;
00108
00109 #endif
00110 goto restart_phi0_with_new_seed;
00111 }
00112
00113 if (mpz_probab_prime_p(x,probab_prime_checks))
00114 {
00115 const unsigned int exponent=mpz_remove(n,n,x);
00116 cout << x << " is factor." << endl;
00117 std::ostringstream comment;
00118 comment << " [" << CRing::name << "0]";
00119 Factorization_to_file << MAL(x,exponent,comment) << flush;
00120 }
00121 else
00122 {
00123 cout << x << " is composite factor." << endl;
00124 mpz_divexact(n,n,x);
00125 if (mpz_probab_prime_p(n,probab_prime_checks))
00126 {
00127 mpz_swap(n,x);
00128 cout << x << " is factor." << endl;
00129 std::ostringstream comment;
00130 comment << " [" << CRing::name << "0/factorswap]";
00131 Factorization_to_file << MAL(x,comment) << flush;
00132
00133 }
00134 else
00135 {
00136 mpz_t saved_n;
00137 mpz_init_set(saved_n,n);
00138 mpz_swap(n,x);
00139 #ifdef VERBOSE_INFO
00140 cout << "Calling " << CRing::name << " (recursive)..." << endl;
00141 #endif
00142 polphi_template<CRing,CRingPhase2>(phase1,phase2);
00143 #ifdef VERBOSE_INFO
00144 cout << "back from recursion. Continuing with " << CRing::name << endl;
00145 #endif
00146 const unsigned int exponent=mpz_remove(saved_n,saved_n,n)+1;
00147 std::ostringstream comment;
00148 comment << " [" << CRing::name << "0]";
00149 if (mpz_probab_prime_p(n,probab_prime_checks))
00150 Factorization_to_file << MAL(n,exponent,comment) << flush;
00151 else
00152 {
00153 comment << " [composite]";
00154 Factorization_to_file << MAL(n,exponent,comment) << flush;
00155 }
00156 mpz_swap(n,saved_n); mpz_clear(saved_n);
00157 }
00158 a.set_startvalue(polphi_seed); i=i_Start;
00159 }
00160 KnownIndices.clear(); PreviousIndex=0;
00161 if (mpz_probab_prime_p(n,probab_prime_checks) || mpz_cmp_ui(n,1)==0) goto done;
00162 }
00163 }
00164 a.pow_mod(2,n);
00165
00166
00167
00168
00169 {
00170
00171 const int D=CRing::SieveIntervalSize;
00172 i=i-(i%D);
00173 cout << CRing::name << "-phase 1" << endl;
00174 while (i<=phase1)
00175 {
00176 #ifdef REACT_ON_SIGUSR
00177 if (USRSignalHandler.got_SIGUSR1()) goto done;
00178 if (USRSignalHandler.got_SIGUSR2()) break;
00179 #endif
00180 b.set(a);
00181
00182 bool sieve[D];
00183 for (int j=1; j<=D; j+=2) sieve[j]=true;
00184 for (int p=5, dp=2; p*p<D+i; p+=dp, dp=6-dp)
00185 for (int j=p-(i%p); j<D; j+=p) sieve[j]=false;
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 for (int j=1, dj=4; j<D; j+=dj, dj=6-dj)
00196 {
00197
00198
00199 if (sieve[j]) a.fast_pow_mod(i+j,n);
00200 }
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214 int old_i=i; i+=D;
00215 #ifdef VERBOSE_NOTICE
00216 cout << "Phase1: " << i << "/" << phase1 << "\r" << flush;
00217 #endif
00218
00219 a.test_gcd(x,n);
00220 if (mpz_cmp_ui(x,1)!=0)
00221 {
00222 cout << "factor found in phase 1: i=" << i << endl;
00223 if (mpz_cmp(x,n)==0)
00224 {
00225 trivialer_factor_found:
00226 cout << "Found factor is trivial, trying to isolate..." << endl;
00227 a.set(b);
00228 i=old_i;
00229 for (int j=1, dj=4; j<D; j+=dj, dj=6-dj) if (sieve[j])
00230 {
00231 a.pow_mod(i+j,n);
00232 a.test_gcd(x,n);
00233 if (mpz_cmp_ui(x,1)!=0)
00234 {
00235 if (mpz_cmp(x,n)==0)
00236 {
00237 cout << "factor is trivial, aborting." << endl;
00238 goto done;
00239 }
00240 else
00241 {
00242 cout << "non-trivial factor isolated at i=" << i+j << endl;
00243 break;
00244 }
00245 }
00246 }
00247 }
00248 if (mpz_probab_prime_p(x,probab_prime_checks))
00249 {
00250 const unsigned int exponent=mpz_remove(n,n,x);
00251 cout << x << " is factor." << endl;
00252 std::ostringstream comment;
00253 comment << " [" << CRing::name << "1]";
00254 Factorization_to_file << MAL(x,exponent,comment) << flush;
00255 }
00256 else
00257 {
00258 mpz_divexact(n,n,x);
00259 cout << x << " is composite factor." << endl;
00260 if (mpz_probab_prime_p(n,probab_prime_checks))
00261 {
00262 cout << n << " is factor." << endl;
00263 std::ostringstream comment;
00264 comment << " [" << CRing::name << "1/factorswap]";
00265 Factorization_to_file << MAL(n,comment) << flush;
00266 mpz_set(n,x); goto trivialer_factor_found;
00267 }
00268 else
00269 {
00270 cout << "Found factor is composite, trying to isolate..." << endl;
00271
00272
00273
00274
00275
00276
00277
00278 mpz_t local_n;
00279 CRing local_a;
00280 mpz_init_set(local_n,x);
00281 local_a.set(b);
00282 for (int j=1, dj=4; j<D; j+=dj, dj=6-dj) if (sieve[j])
00283 {
00284 local_a.pow_mod(old_i+j,local_n);
00285 local_a.test_gcd(x,local_n);
00286 if (mpz_cmp_ui(x,1)!=0)
00287 {
00288 mpz_divexact(local_n,local_n,x);
00289 const unsigned int exponent=mpz_remove(n,n,x)+1;
00290 cout << x << " is factor." << endl;
00291 if (mpz_probab_prime_p(x,probab_prime_checks))
00292 {
00293 cout << "factor isolated at i=" << old_i+j << endl;
00294 std::ostringstream comment;
00295 comment << " [" << CRing::name << "1]";
00296 Factorization_to_file << MAL(x,exponent,comment) << flush;
00297 }
00298 else
00299 {
00300 cout << "factor still composite at i=" << old_i+j << endl;
00301 std::ostringstream comment;
00302 comment << " [" << CRing::name << "1] [composite]";
00303 Factorization_to_file << MAL(x,exponent,comment) << flush;
00304 }
00305 if (mpz_cmp_ui(local_n,1)==0) break;
00306 }
00307 }
00308 if (mpz_cmp_ui(local_n,1)!=0)
00309 {
00310 MARK; cerr << "Error: local_n should be 1 now!" << endl; exit(1);
00311 }
00312 mpz_clear(local_n);
00313 }
00314 }
00315 if (mpz_probab_prime_p(n,probab_prime_checks) || mpz_cmp_ui(n,1)==0) goto done;
00316 }
00317 }
00318 i--; d=4; while (!is_prime(i)) { i-=d; d=6-d; }
00319 #ifdef VERBOSE_NOTICE
00320 cout << CRing::name << "-phase 1, up to " << i << endl;
00321 #endif
00322 if (i>=phase2) goto done;
00323 }
00324
00325
00326
00327
00328
00329 #if CONTINUATION_METHOD > 0
00330
00331 {
00332 double ii=i;
00333
00334
00335
00336 cout << CRing::name << "-phase 2 (improved with pairing & fft)" << endl;
00337
00338
00339 unsigned int pg=256, pms=2*1050;
00340 get_fft_parameter(pg,pms,phase2,n);
00341
00342 const unsigned int D = pms;
00343 const unsigned int Polynom_Index = pg;
00344
00345
00346
00347
00348 const polynomial::TPolynom collector = new mpz_t[Polynom_Index];
00349 for (unsigned int j=0; j<Polynom_Index; ++j) mpz_init(collector[j]);
00350
00351
00352 CRingPhase2 a_phase2(a);
00353
00354
00355 unsigned int Polynom_Index_i = 0;
00356 for (unsigned int j=1; j<D/2; j+=2)
00357 {
00358 if (gcd(j,D)==1)
00359 {
00360 a_phase2.get_polynomdef_point(x);
00361 mpz_set(collector[Polynom_Index_i++],x);
00362 }
00363 a_phase2.calc_polynomdef_next_point();
00364 }
00365
00366 #ifdef VERBOSE_INFO
00367 cout << "original size of polynomial: " << Polynom_Index_i << endl;
00368 #endif
00369 polynomial::TPolynom Polynom = NULL;
00370
00371
00372 {
00373 unsigned int anz_Punkte=Polynom_Index_i;
00374 while (anz_Punkte<Polynom_Index) mpz_set_ui(collector[anz_Punkte++],0);
00375 Polynom_Index_i=polynomial::construct_polynomial_from_roots(Polynom,collector,anz_Punkte,n);
00376 #ifdef VERBOSE_INFO
00377 cout << "polynomial created. size of polynomial = " << Polynom_Index_i << endl;
00378 #endif
00379 }
00380
00381 ii=floor(ii/D)*D;
00382 a_phase2.calc_EvalStartingPoint(D,ii);
00383 ii-=D/2;
00384 while (ii<=phase2)
00385 {
00386 #ifdef VERBOSE_NOTICE
00387 cout << "Phase2: " << std::setprecision(0) << ii << "/" << phase2 << endl << flush;
00388 #endif
00389 #ifdef VERBOSE_INFO
00390 cout << "Computing and collecting values..." << endl;
00391 #endif
00392
00393
00394 unsigned int collector_index=0;
00395 while (collector_index<(Polynom_Index_i-1)/2)
00396 {
00397 a_phase2.get_point_and_calc_next_point(x);
00398 mpz_set(collector[collector_index++],x);
00399
00400 ii+=D;
00401 }
00402
00403 #ifdef REACT_ON_SIGUSR
00404 if (USRSignalHandler.got_SIGUSR1()) goto done_phase2;
00405 if (USRSignalHandler.got_SIGUSR2()) continue;
00406 #endif
00407
00408
00409 #ifdef VERBOSE_NOTICE
00410 cout << "Phase2: " << std::setprecision(0) << ii << "/" << phase2 << endl << flush;
00411 #endif
00412 #ifdef VERBOSE_INFO
00413 cout << "starting multipoint polynomial evaluation" << endl;
00414 cout << "Parameter: polynomial size: " << Polynom_Index_i
00415 << ", points to evaluate: " << collector_index << endl;
00416 #endif
00417 polynomial::multipoint_eval(collector,Polynom,Polynom_Index_i,collector,collector_index,n);
00418 #ifdef VERBOSE_INFO
00419 cout << "polynomial evaluation finished, computing gcd." << endl;
00420 #endif
00421 mpz_set(y,collector[0]);
00422 for (unsigned int j=1; j<collector_index; ++j)
00423 { mpz_mul(y,y,collector[j]); mpz_mod(y,y,n); }
00424
00425 mpz_gcd(x,y,n);
00426 if (mpz_cmp_ui(x,1)==0) continue;
00427
00428
00429
00430
00431
00432
00433
00434
00435 for (unsigned int j=0; j<collector_index; ++j)
00436 {
00437 mpz_gcd(x,collector[j],n);
00438 if (mpz_cmp_ui(x,1)!=0)
00439 {
00440 cout << "factor found in phase 2: i=" << ii << endl;
00441 if (mpz_probab_prime_p(x,probab_prime_checks))
00442 {
00443 const unsigned int exponent=mpz_remove(n,n,x);
00444 cout << x << " is factor." << endl;
00445 std::ostringstream comment;
00446 comment << " [" << CRing::name << "2]";
00447 Factorization_to_file << MAL(x,exponent,comment) << flush;
00448 }
00449 else
00450 {
00451 mpz_divexact(n,n,x);
00452 cout << x << " is composite factor." << endl;
00453 if (mpz_probab_prime_p(n,probab_prime_checks))
00454 {
00455 mpz_swap(n,x);
00456 cout << x << " is factor." << endl;
00457 std::ostringstream comment;
00458 comment << " [" << CRing::name << "2/factorswap]";
00459 Factorization_to_file << MAL(x,comment) << flush;
00460 }
00461 else
00462 {
00463 if (mpz_cmp_ui(n,1)==0)
00464 {
00465
00466 mpz_swap(n,x);
00467 mpz_set_d(y,ii);
00468 Factorization_to_file << "! [" << CRing::name << "2/trivial! i=" << y << "]" << flush;
00469 goto done_phase2;
00470 }
00471 else
00472 {
00473 std::ostringstream comment;
00474 comment << " [" << CRing::name << "2] [composite]";
00475 Factorization_to_file << MAL(x,comment) << flush;
00476 }
00477 }
00478 }
00479 if (mpz_probab_prime_p(n,probab_prime_checks)) goto done_phase2;
00480 if (mpz_cmp_ui(n,1)==0) goto done_phase2;
00481
00482 for (unsigned int j=0; j<Polynom_Index_i; ++j) mpz_mod(Polynom[j],Polynom[j],n);
00483 }
00484 }
00485
00486 }
00487
00488 done_phase2:
00489 cout << CRing::name << "-phase 2, up to " << std::setprecision(0) << ii << endl;
00490 for (unsigned int j=0; j<Polynom_Index; ++j) mpz_clear(collector[j]);
00491 delete [] collector;
00492 for (unsigned int j=0; j<Polynom_Index_i; ++j) mpz_clear(Polynom[j]);
00493 delete [] Polynom;
00494 }
00495 #else
00496 #warning continuation disabled, because CONTINUATION_METHOD 0 not implemented
00497 #endif
00498
00499 done:
00500 mpz_clear(x); mpz_clear(y);
00501 #ifdef VERBOSE_INFO
00502 cout << CRing::name << " method finished..." << endl;
00503 #endif
00504 }
00505
00506 #endif