mpqsMultiplier.cc

Go to the documentation of this file.

#include "utils.H"
#include "mpz_wrapper.H"
#include "modulo.H"
#include <cmath>

// switch for method to compute a multiplier
#define BEWERTUNG_MIT_POTENZEN
//#define BEWERTUNG_MIT_DIRTY_SIEVING


void determine_best_MPQS_Multiplier(const mpz_t n, mpz_t kN, int &new_MPQS_Multiplier)
{
  // Now we're going to determine a convenient multiplier.
  // result: new_MPQS_Multiplier and kN = new_MPQS_Multiplier*n

  /* Remark: In this implementation for the multiple polynomial quadratic sieve
     the constraint kN=1 (mod 4) must be met, too.
     It is also advisable to fulfil kN=1 (mod 8) to assure that 2 is member of the factorbase.
     (cf. Silverman paper) */

  using numtheory::is_prime;
  using namespace my_mpz_wrapper;

  // determine the multiplier:
  const int P=2; // biggest prime, up to which members of the factorbase are consecutive primes.
                 // we need at least P=2.
  int bound_k = 1000000; /* maximal die ersten bis zum Vorfaktor bound_k den Test
                               durchführen. (Während des Tests wird diese Schranke dann
                               höchstens noch verkleinert.) */
  const int FB_SIM_Size = 100; /* Faktorbasisgröße, mit der Simulation durchgeführt werden soll. */
  
  // 1. Es sollen alle Primzahlen von 2 bis P in der statischen Faktorbasis enthalten sein.
  //    Dafür wird in Kauf genommen, dass der Vorfaktor etwas größer wird...
  // 2. Maximal bis "bound_k" wird für Vorfaktoren simuliert, wie gut diese Vorfaktoren abschneiden.
  //    Hierfür werden Teilbarkeitswahrscheinlichkeiten verwendet und die Größenordnung
  //    der zu siebenden Zahl anhand des Vorfaktors berücksichtigt.
  // 3. Wird ein neuer besserer Vorfaktor entdeckt, dann wird "bound_k" neu kalkuliert.
  // 4. Der hierbei am besten abschneidende Vorfaktor wird schließlich verwendet.

#ifdef VERBOSE_INFO  
  cout << "Determining a multiplier for MPQS" << endl;
#endif

  if (mpz_sizeinbase(n,10)<50)
   bound_k=100; // für kleine Zahlen begrenzen, da sonst die Vorfaktorermittlung länger als Faktorisierung dauert...
  
  // zunächst die theoretisch (erreichbare) optimale Faktorbasis als Berechnungsgrundlage
  // für Optimalwert ermitteln:
  double optimal_w, good_w;
  int p; // für Primzahltest
  
  p=2; optimal_w= -2*log(2);
  while (p<P)
    { 
      while (!is_prime(p)) p++;
#ifdef BEWERTUNG_MIT_POTENZEN
      optimal_w -= 2*log(p)/(p-1);
#else
      optimal_w -= 2*log(p)/p;
#endif
      p++;
    }
  // bis hierher alle Primzahlen in Folge in simulierter FB
  // für theoretisch optimale FB müssen sämtliche Primzahlen in Folge in FB liegen
  for (int z=0; z<FB_SIM_Size; p++)
    { 
      while (!is_prime(p)) p++;
      z++;
#ifdef BEWERTUNG_MIT_POTENZEN
      optimal_w -= 2*log(p)/(p-1);
#else
      optimal_w -= 2*log(p)/p;
#endif
    }

  // Nun noch eine realistischere Bewertungsgrundlage:
  good_w=0.85*optimal_w;

#ifdef VERBOSE_INFO  
  cout << "theoretically optimal valuation: " << std::setprecision(10) << optimal_w << endl;
  cout << "theoretically good valuation   : " << good_w << endl; 
  cout << "Starting simulation of MPQS multipliers..." << endl;
#endif
  
  // ab hier die "wirkliche" Vorfaktorermittlung

  // Besonderheit: Dirty-Sieving wird korrekt simuliert, d.h.
  // für die ersten "SieveControl::FBLowestStartIndex" Primzahlen wird
  // die Trefferwahrscheinlichkeit um 1/p vermindert.  


  // Feld für sortierte Aufnahme der "sort_bis"+1 besten Vorfaktoren
  struct tFeldeintrag
   {
     double w; int k;
     inline tFeldeintrag() : w(0.0), k(0) { }
     inline ~tFeldeintrag() { }
   };
  const short int sort_bis = 5;    
  tFeldeintrag Feld[sort_bis]; /* Feld wird automatisch mit Nullen initialisiert
                                  -> für Bewertung werden daher grundsätzlich
                                  negative Werte vorausgesetzt */

  double best_w = 1e20;
  double akt_w;
  int best_k = 1;
  int k=1; // testweiser Vorfaktor
  int z=0; // Index der Primzahl in der Faktorbasis
  mpz_t x; mpz_init(x);
  
  while (k<=bound_k)
    {
      mpz_mul_ui(kN,n,k); // initialize kN
      
      // is 2 member of the factorbase?
      if (mpz_remainder_ui(kN,8)!=1) goto next;

      akt_w=-2*log(2); // Initialisierung der Bewertung für neuen Vorfaktor
      
      p=2; z=1;
      while (p<P)
        {
          do p++; while(!is_prime(p));
          mpz_set_ui(x,p); if (mpz_legendre(kN,x)==-1) goto next; // goto next -> reject this multiplier
          z++;

          // Aktualisierung der Bewertung:
#ifdef BEWERTUNG_MIT_POTENZEN
#ifdef BEWERTUNG_MIT_DIRTY_SIEVING
          if (z<SieveControl::FBLowestStartIndex)
            akt_w -= (k%p!=0)?  2*log(p)/(p-1) - 1.75*log(p)/p: 0.1*log(p)/p;
          else
            akt_w -= (k%p!=0)?  2*log(p)/(p-1) : 1*log(p)/p;
#else // ohne BEWERTUNG_MIT_DIRTY_SIEVING
          akt_w -= (k%p!=0)?  2*log(p)/(p-1) : 1*log(p)/p;        
#endif
#else
#ifdef BEWERTUNG_MIT_DIRTY_SIEVING
          if (z<SieveControl::FBLowestStartIndex)
            akt_w -= 0.5*((k%p!=0)?  2*log(p)/p : log(p)/p);
          else
            akt_w -= (k%p!=0)?  2*log(p)/p : 1*log(p)/p;
#else // ohne BEWERTUNG_MIT_DIRTY_SIEVING
          akt_w -= (k%p!=0)?  2*log(p)/p : 1*log(p)/p;
#endif
#endif
          
        }
      
      // Vorfaktor prinzipiell geeignet, nun noch weiter bewerten:
      {
        for ( ; z<FB_SIM_Size; z++)
          {
            do
              {
                do p++; while(!is_prime(p)); /* nächste Primzahl */
                mpz_set_ui(x,p);
              }
            while (mpz_legendre(kN,x)==-1); // Primzahl gültig für Primzahlbasis?


            // Aktualisierung der Bewertung:
#ifdef BEWERTUNG_MIT_POTENZEN
#ifdef BEWERTUNG_MIT_DIRTY_SIEVING
          if (z<SieveControl::FBLowestStartIndex)
            akt_w -= (k%p!=0)?  2*log(p)/(p-1) - 1.75*log(p)/p : 0.1*log(p)/p;
          else
            akt_w -= (k%p!=0)?  2*log(p)/(p-1) : 1*log(p)/p;
#else // ohne BEWERTUNG_MIT_DIRTY_SIEVING
          akt_w -= (k%p!=0)?  2*log(p)/(p-1) : 1*log(p)/p;        
#endif
#else
#ifdef BEWERTUNG_MIT_DIRTY_SIEVING
          if (z<SieveControl::FBLowestStartIndex)
            akt_w -= 0.5*((k%p!=0)?  2*log(p)/p : log(p)/p);
          else
            akt_w -= (k%p!=0)?  2*log(p)/p : 1*log(p)/p;
#else // ohne BEWERTUNG_MIT_DIRTY_SIEVING
          akt_w -= (k%p!=0)?  2*log(p)/p : 1*log(p)/p;
#endif
#endif
            
          }

        akt_w += 0.5*log(k); // add threshold

        // nun noch den Vorfaktor mit Bewertung ausgeben
        //cout << k << " mit Bewertung " << akt_w << endl;

        { 
          // for statistical reasons, collect the best five multipliers:
          signed short int i=sort_bis-1;
          if (akt_w<Feld[i].w)
            {
              while (i>0 && akt_w<Feld[i-1].w) { Feld[i]=Feld[i-1]; --i; }
              Feld[i].w=akt_w; Feld[i].k=k;
            }
        }
        
        if (akt_w<best_w)
          {
#ifdef VERBOSE
            cout << "best multiplier so far: " <<  k 
                 << " valuation: " << akt_w << endl;
#endif
            best_w=akt_w; best_k=k;
            int computed_bound_k;
            if (akt_w>good_w)
                 computed_bound_k = static_cast<int>(ceil(exp(2*(akt_w-good_w))));
            else computed_bound_k = 0;
            if (computed_bound_k<bound_k) bound_k=computed_bound_k;
#ifdef VERBOSE
            cout << "New (practical) bound: " << bound_k
                 << " (theoretical bound: " << computed_bound_k << ")" <<  endl;
#endif
          }
      }
      
      
    next: // next multiplier
#if 0
      do k++; while(!is_prime(k)); // in case that only prime numbers are allowed as multipliers
#else
      // alternative: allow more multipliers (constraint: squarefree, positive integers)
      k++;
      { // is k squarefree?
        int square=4, i=5;
        while (square<=k)
          {
            if (k%square==0) goto next; // contains a square
            square+=i; i+=2;
          }
        // okay, k is squarefree...
      }
#endif
    }

#ifdef VERBOSE_INFO
  cout << endl;
  cout << "The " << sort_bis << " best MPQS multipliers are" << endl;
  for (short int i=0; i<sort_bis; i++)
    cout << i+1 << ". " << Feld[i].k << " with valuation " << Feld[i].w << endl;
  cout << endl;
#endif
 
#if 0
  // If this block is activated, one can type in a multiplier manually.
  cout << "Please input a multiplier: ";
  cin >> best_k;
#endif

  new_MPQS_Multiplier=best_k; // this is the new MPQS_Multiplier
#ifdef VERBOSE_NOTICE
  cout << "MPQS multiplier set to " << best_k << endl;
#endif
  mpz_mul_ui(kN,n,best_k); // compute best kN
  mpz_clear(x);
}

---