#ifndef polphi_template_HEADER
#define polphi_template_HEADER
/*! @file
* @brief
* template definition for pollard-phi-like factoring algorithms
*/
#include <vector>
#include <sstream>
extern void get_fft_parameter(unsigned int &pg_returnvalue, unsigned int &pms_returnvalue,
const double phase2, const mpz_t N);
template <class CRing, class CRingPhase2>
void polphi_template(const int phase1, const double phase2)
{
// Pollard-Phi method modified in three phases:
// Phase 0: small numbers up to a constant upper limit
// Phase 1: prime numbers up to "phase1" as upper limit
// [ factor will be detected, if a prime factor P minus 1 is a product of some of these prime numbers ]
// Phase 2: prime numbers up to "phase2" as upper limit
// [ factor will be detected, if the only one missing prime factor of P-1 in phase 1 is within phase 2 ]
// (improved standard continuation with pairing & fft)
using numtheory::is_prime;
using numtheory::gcd;
mpz_t x,y;
mpz_init(x); mpz_init(y);
CRing a,b;
std::vector<int> KnownIndices;
int PreviousIndex = 0;
unsigned int polphi_seed = 0;
const int i_Start = 2;
int i;
short int d;
restart_phi0_with_new_seed:
KnownIndices.clear(); PreviousIndex=0;
++polphi_seed;
//
// Phase 0: small numbers (and powers of small prime numbers)
//
cout << "--------------------------------------------------" << endl;
cout << CRing::name << "-phase 0" << endl;
if (!a.set_startvalue(polphi_seed)) goto done; // seed wasn't accepted
for (i=i_Start; i<10000; i++)
{
a.pow_mod(i,n);
retry:
a.test_gcd(x,n);
if (mpz_cmp_ui(x,1)!=0) // factor found?
{
cout << "factor found in phase 0: i=" << i << endl;
if (mpz_cmp(x,n)==0) // gefundener Faktor = n ?
{
cout << "factor is trivial." << endl;
int h=i; int j=2;
while (h>j) // determine biggest prime factor
{
while(h>j && h%j==0) h/=j;
++j;
}
if (KnownIndices.empty()) KnownIndices.push_back(h);
else // was it added in previous run?
if (PreviousIndex!=i || KnownIndices.back()!=h) KnownIndices.push_back(h);
else
{
// then -- should the occasion arise -- remove it!
if (KnownIndices.back()!=i/h) KnownIndices.pop_back();
// and replace it by another divisor
KnownIndices.push_back(i/h);
}
PreviousIndex=i;
#ifdef VERBOSE_INFO
cout << "KnownIndices:";
for (std::vector<int>::const_iterator p=KnownIndices.begin(); p!=KnownIndices.end(); ++p) cout << " " << *p;
cout << endl;
#endif
if (KnownIndices.empty() || KnownIndices.back()==1)
{
#ifdef VERBOSE_INFO
MARK; cout << "Giving up for " << n << endl;
//{ char ch; cin >> ch; }
#endif
goto restart_phi0_with_new_seed; // trivial factor -> try to roll up factors differently or abort this procedure
}
if (KnownIndices.size()<100)
{
a.set_startvalue(polphi_seed);
for (std::vector<int>::const_iterator p=KnownIndices.begin(); p!=KnownIndices.end(); ++p) a.pow_mod(*p,n);
cout << "Trying again..." << endl;
i=i_Start; goto retry;
}
#ifdef VERBOSE_WARN
MARK; cout << "Aborting for " << n << endl;
//{ char ch; cin >> ch; }
#endif
goto restart_phi0_with_new_seed; // trivial factor -> try to roll up factors differently or abort this procedure
}
if (mpz_probab_prime_p(x,probab_prime_checks))
{
const unsigned int exponent=mpz_remove(n,n,x);
cout << x << " is factor." << endl;
std::ostringstream comment;
comment << " [" << CRing::name << "0]";
Factorization_to_file << MAL(x,exponent,comment) << flush;
}
else
{
cout << x << " is composite factor." << endl;
mpz_divexact(n,n,x);
if (mpz_probab_prime_p(n,probab_prime_checks))
{
mpz_swap(n,x);
cout << x << " is factor." << endl;
std::ostringstream comment;
comment << " [" << CRing::name << "0/factorswap]";
Factorization_to_file << MAL(x,comment) << flush;
//goto done;
}
else
{
mpz_t saved_n;
mpz_init_set(saved_n,n);
mpz_swap(n,x);
#ifdef VERBOSE_INFO
cout << "Calling " << CRing::name << " (recursive)..." << endl;
#endif
polphi_template<CRing,CRingPhase2>(phase1,phase2);
#ifdef VERBOSE_INFO
cout << "back from recursion. Continuing with " << CRing::name << endl;
#endif
const unsigned int exponent=mpz_remove(saved_n,saved_n,n)+1;
std::ostringstream comment;
comment << " [" << CRing::name << "0]";
if (mpz_probab_prime_p(n,probab_prime_checks))
Factorization_to_file << MAL(n,exponent,comment) << flush;
else
{
comment << " [composite]";
Factorization_to_file << MAL(n,exponent,comment) << flush;
}
mpz_swap(n,saved_n); mpz_clear(saved_n);
}
a.set_startvalue(polphi_seed); i=i_Start;
}
KnownIndices.clear(); PreviousIndex=0;
if (mpz_probab_prime_p(n,probab_prime_checks) || mpz_cmp_ui(n,1)==0) goto done;
}
}
a.pow_mod(2,n); // guarantee that index is an even value
//
// Phase 1: prime numbers (without powers)
//
{
// Phase 1: prime numbers (without powers)
const int D=CRing::SieveIntervalSize;
i=i-(i%D); // i with constraint i=5 (mod 6) for fast calculation of prime numbers
cout << CRing::name << "-phase 1" << endl;
while (i<=phase1)
{
#ifdef REACT_ON_SIGUSR
if (USRSignalHandler.got_SIGUSR1()) goto done;
if (USRSignalHandler.got_SIGUSR2()) break;
#endif
b.set(a); // memorize old value (to be able to roll up, if necessary)
// sieving prime numbers
bool sieve[D];
for (int j=1; j<=D; j+=2) sieve[j]=true; // initialize sieve
for (int p=5, dp=2; p*p<D+i; p+=dp, dp=6-dp)
for (int j=p-(i%p); j<D; j+=p) sieve[j]=false;
/*
Note:
a) p*p<D+i is a little bit too strong as constraint for
detecting prime numbers, but that doesn't matter because of Phase 0.
(-> Prime numbers below sqrt(D) are not detected)
b) p-(i%p) is p, if p divides i; (this is unwanted!)
but: we don't need Index 0 later on, anyway. So it doesn't matter.
*/
for (int j=1, dj=4; j<D; j+=dj, dj=6-dj)
{
//if (sieve[j]!=is_prime(i+j))
// { cerr << "inconsistency for: " << i+j << endl; }
if (sieve[j]) a.fast_pow_mod(i+j,n);
}
/* Bemerkung:
Nicht durch die Literatur verwirren lassen:
Das akkumulierte Produkt der (a-1)-Werte muß natürlich
nicht gebildet werden. Das wäre eine total überflüssige
Aktion, da ein Teiler auch bei nachfolgenden a's erhalten
bleibt (aufgrund des kleinen Fermatschen Satzes!).
Falls sich andererseits mehrere Teiler akkumulieren, so kann
auch das akkumulierte Produkt dies nicht verhindern...
*/
// the loop avoids frequent gcd-computations;
// in seldom cases multiple factors may be found, which cannot be separated.
// -> Phi-Method fails in these cases (if necessary: reduce loop size D)
int old_i=i; i+=D;
#ifdef VERBOSE_NOTICE
cout << "Phase1: " << i << "/" << phase1 << "\r" << flush;
#endif
a.test_gcd(x,n);
if (mpz_cmp_ui(x,1)!=0)
{
cout << "factor found in phase 1: i=" << i << endl;
if (mpz_cmp(x,n)==0) // is found factor = n ?
{
trivialer_factor_found:
cout << "Found factor is trivial, trying to isolate..." << endl;
a.set(b); // Basis vom Intervallanfang
i=old_i; // zurück zum Intervallanfang
for (int j=1, dj=4; j<D; j+=dj, dj=6-dj) if (sieve[j])
{
a.pow_mod(i+j,n);
a.test_gcd(x,n);
if (mpz_cmp_ui(x,1)!=0)
{
if (mpz_cmp(x,n)==0) // is found factor = n ?
{
cout << "factor is trivial, aborting." << endl;
goto done; // trivial factor -> try to roll up in another fashion or abort this procedure
}
else
{
cout << "non-trivial factor isolated at i=" << i+j << endl;
break; // factor isolated...
}
}
}
}
if (mpz_probab_prime_p(x,probab_prime_checks))
{
const unsigned int exponent=mpz_remove(n,n,x);
cout << x << " is factor." << endl;
std::ostringstream comment;
comment << " [" << CRing::name << "1]";
Factorization_to_file << MAL(x,exponent,comment) << flush;
}
else
{
mpz_divexact(n,n,x);
cout << x << " is composite factor." << endl;
if (mpz_probab_prime_p(n,probab_prime_checks))
{
cout << n << " is factor." << endl;
std::ostringstream comment;
comment << " [" << CRing::name << "1/factorswap]";
Factorization_to_file << MAL(n,comment) << flush;
mpz_set(n,x); goto trivialer_factor_found;
}
else
{
cout << "Found factor is composite, trying to isolate..." << endl;
// The found factor is composite and has already been divided out of our number n.
// We could multiply him back again, but we can do much better:
// As we are only interested in splitting our composite factor, we can
// repeat the search in the previous interval, but this time modulo the composite factor!
// For this reason we need local variables for the base "a" and our "n", because
// their original values are needed in the main loop, so we shouldn't clobber them.
mpz_t local_n;
CRing local_a;
mpz_init_set(local_n,x);
local_a.set(b); // Basis vom Intervallanfang
for (int j=1, dj=4; j<D; j+=dj, dj=6-dj) if (sieve[j])
{
local_a.pow_mod(old_i+j,local_n);
local_a.test_gcd(x,local_n);
if (mpz_cmp_ui(x,1)!=0)
{
mpz_divexact(local_n,local_n,x);
const unsigned int exponent=mpz_remove(n,n,x)+1;
cout << x << " is factor." << endl;
if (mpz_probab_prime_p(x,probab_prime_checks))
{
cout << "factor isolated at i=" << old_i+j << endl;
std::ostringstream comment;
comment << " [" << CRing::name << "1]";
Factorization_to_file << MAL(x,exponent,comment) << flush;
}
else
{
cout << "factor still composite at i=" << old_i+j << endl;
std::ostringstream comment;
comment << " [" << CRing::name << "1] [composite]";
Factorization_to_file << MAL(x,exponent,comment) << flush;
}
if (mpz_cmp_ui(local_n,1)==0) break;
}
}
if (mpz_cmp_ui(local_n,1)!=0)
{
MARK; cerr << "Error: local_n should be 1 now!" << endl; exit(1);
}
mpz_clear(local_n);
}
}
if (mpz_probab_prime_p(n,probab_prime_checks) || mpz_cmp_ui(n,1)==0) goto done;
}
}
i--; d=4; while (!is_prime(i)) { i-=d; d=6-d; } // adjust "i" to previous prime number
#ifdef VERBOSE_NOTICE
cout << CRing::name << "-phase 1, up to " << i << endl;
#endif
if (i>=phase2) goto done;
}
//
// Phase 2
//
#if CONTINUATION_METHOD > 0
/* improved standard continuation with pairing & fft */
{
double ii=i;
//const double i = -1; // uncomment this just to be sure that i isn't used anymore...
// Phase 2: isolated (prime) number: p-1 has to be completely decomposable before, except for the one factor to detect within this phase
cout << CRing::name << "-phase 2 (improved with pairing & fft)" << endl;
// for info about these values: see elliptic curves, fft-continuation!
unsigned int pg=256, pms=2*1050;
get_fft_parameter(pg,pms,phase2,n); // get parameter (pg,pms) for "phase2"
const unsigned int D = pms;
const unsigned int Polynom_Index = pg;
//cout << "size of polynomial: " << Polynom_Index << endl;
//cout << "single interval size:" << D << endl;
const polynomial::TPolynom collector = new mpz_t[Polynom_Index];
for (unsigned int j=0; j<Polynom_Index; ++j) mpz_init(collector[j]);
// "a" is the initial value for Phase 2
CRingPhase2 a_phase2(a);
// Precomputing most common used values:
unsigned int Polynom_Index_i = 0;
for (unsigned int j=1; j<D/2; j+=2)
{
if (gcd(j,D)==1)
{
a_phase2.get_polynomdef_point(x);
mpz_set(collector[Polynom_Index_i++],x);
}
a_phase2.calc_polynomdef_next_point();
}
#ifdef VERBOSE_INFO
cout << "original size of polynomial: " << Polynom_Index_i << endl;
#endif
polynomial::TPolynom Polynom = NULL; // polynomial we're going to construct
// fill up with zeros up to the next power of two
// Note: It would be better to add other points outside our search interval.
{
unsigned int anz_Punkte=Polynom_Index_i;
while (anz_Punkte<Polynom_Index) mpz_set_ui(collector[anz_Punkte++],0);
Polynom_Index_i=polynomial::construct_polynomial_from_roots(Polynom,collector,anz_Punkte,n);
#ifdef VERBOSE_INFO
cout << "polynomial created. size of polynomial = " << Polynom_Index_i << endl;
#endif
}
ii=floor(ii/D)*D; // that is: ii=ii-(ii%D); // determine starting point
a_phase2.calc_EvalStartingPoint(D,ii);
ii-=D/2; // normalize interval from [i-D/2,i+D/2] to [i,i+D] (to ease notation)
while (ii<=phase2)
{
#ifdef VERBOSE_NOTICE
cout << "Phase2: " << std::setprecision(0) << ii << "/" << phase2 << endl << flush;
#endif
#ifdef VERBOSE_INFO
cout << "Computing and collecting values..." << endl;
#endif
// collecting values
unsigned int collector_index=0;
while (collector_index<(Polynom_Index_i-1)/2)
{
a_phase2.get_point_and_calc_next_point(x);
mpz_set(collector[collector_index++],x);
// compute values for next interval...
ii+=D; // new interval (up to "i" the interval has been calculated)
}
#ifdef REACT_ON_SIGUSR
if (USRSignalHandler.got_SIGUSR1()) goto done_phase2;
if (USRSignalHandler.got_SIGUSR2()) continue;
#endif
// start polynomial evaluation
#ifdef VERBOSE_NOTICE
cout << "Phase2: " << std::setprecision(0) << ii << "/" << phase2 << endl << flush;
#endif
#ifdef VERBOSE_INFO
cout << "starting multipoint polynomial evaluation" << endl;
cout << "Parameter: polynomial size: " << Polynom_Index_i
<< ", points to evaluate: " << collector_index << endl;
#endif
polynomial::multipoint_eval(collector,Polynom,Polynom_Index_i,collector,collector_index,n);
#ifdef VERBOSE_INFO
cout << "polynomial evaluation finished, computing gcd." << endl;
#endif
mpz_set(y,collector[0]);
for (unsigned int j=1; j<collector_index; ++j)
{ mpz_mul(y,y,collector[j]); mpz_mod(y,y,n); }
// nun (erst) den ggT ermitteln
mpz_gcd(x,y,n);
if (mpz_cmp_ui(x,1)==0) continue;
// now one or more factor have been found
// possible problem: bundled factors which have to be separated
// Da es aber nur log-viele Teiler geben kann, ist dieser
// Abschnitt nicht zeitkritisch; also stört es wenig, alle
// Werte nochmal durchzugehen (dies senkt die Wahrscheinlichkeit,
// akkumulierte Teiler zu erhalten).
for (unsigned int j=0; j<collector_index; ++j)
{
mpz_gcd(x,collector[j],n);
if (mpz_cmp_ui(x,1)!=0)
{
cout << "factor found in phase 2: i=" << ii << endl;
if (mpz_probab_prime_p(x,probab_prime_checks))
{
const unsigned int exponent=mpz_remove(n,n,x);
cout << x << " is factor." << endl;
std::ostringstream comment;
comment << " [" << CRing::name << "2]";
Factorization_to_file << MAL(x,exponent,comment) << flush;
}
else
{
mpz_divexact(n,n,x);
cout << x << " is composite factor." << endl;
if (mpz_probab_prime_p(n,probab_prime_checks))
{
mpz_swap(n,x);
cout << x << " is factor." << endl;
std::ostringstream comment;
comment << " [" << CRing::name << "2/factorswap]";
Factorization_to_file << MAL(x,comment) << flush;
}
else
{
if (mpz_cmp_ui(n,1)==0)
{
// this is really interesting!
mpz_swap(n,x);
mpz_set_d(y,ii); // short hack to get a well-formatted ii
Factorization_to_file << "! [" << CRing::name << "2/trivial! i=" << y << "]" << flush;
goto done_phase2;
}
else
{
std::ostringstream comment;
comment << " [" << CRing::name << "2] [composite]";
Factorization_to_file << MAL(x,comment) << flush;
}
}
}
if (mpz_probab_prime_p(n,probab_prime_checks)) goto done_phase2;
if (mpz_cmp_ui(n,1)==0) goto done_phase2;
// reduce polynomial modulo n:
for (unsigned int j=0; j<Polynom_Index_i; ++j) mpz_mod(Polynom[j],Polynom[j],n);
}
}
}
done_phase2:
cout << CRing::name << "-phase 2, up to " << std::setprecision(0) << ii << endl;
for (unsigned int j=0; j<Polynom_Index; ++j) mpz_clear(collector[j]);
delete [] collector;
for (unsigned int j=0; j<Polynom_Index_i; ++j) mpz_clear(Polynom[j]);
delete [] Polynom;
}
#else
#warning continuation disabled, because CONTINUATION_METHOD 0 not implemented
#endif
done:
mpz_clear(x); mpz_clear(y);
#ifdef VERBOSE_INFO
cout << CRing::name << " method finished..." << endl;
#endif
}
#endif