polphi_template.H

Go to the documentation of this file.

#ifndef polphi_template_HEADER
#define polphi_template_HEADER

#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

---