qsieve.H

Go to the documentation of this file.

#ifndef qsieve_header
#define qsieve_header

#include "qsieve-fwd.H"
#include "utils.H"

#ifdef CYGWIN_COMPAT
 extern "C"
 {
   #include <malloc.h>
 }
#endif


// use wrapped C++ stream I/O for mpz_t (multiple precision numbers)
#include "mpz_wrapper.H"
using namespace my_mpz_wrapper;

#include "Tfactor.H"
#include "TinyVector.H"
typedef CTinyVector<StaticFactorbaseSettings::FBsizetype> CTinyFBsizetypeVector;

#include "myBitString.H"

using std::streampos;


class CStreamEncoder
{
 private:
  std::ostream &out;
  int prev_value; // previous value to write (but not yet written because of incomplete encoding)
 public:
  CStreamEncoder(std::ostream &out_) : out(out_), prev_value(0) { }

  void PutValue(const int w);
};

class CStreamDecoder
{
 private:
  std::istream &in;
  int ahead_value; // already decoded value (will be returned on next read)
 public:
  CStreamDecoder(std::istream &in_) : in(in_), ahead_value(0) { }

  int GetValue();
};

class CProvideHelperVariables : private ForbidAssignment, protected StaticFactorbaseSettings
{
 // it's not nice, but better than providing temporary mpz_t x,y
 // in the global namespace: threadsafe & restricted to the class
 // And we have no need to change the source code at 1000 places...
 protected:
  mpz_t x,y; // two helper variables
 public:
  CProvideHelperVariables() { mpz_init(x); mpz_init(y); }
  ~CProvideHelperVariables() { mpz_clear(x); mpz_clear(y); }
};

class CRelation : protected CProvideHelperVariables
/* CRelation manages a relation of the form
       Delta^2 = product of factors (mod N)
   Our aim is to reduce the set of factors by combining relations.
   If two factors are equal, their product will be a square and can
   thus be eliminated (by taking the factor into Delta).

   If all factors are eliminated, the product on the right side will evaluate
   to 1. -> The system of equations is solved.

   Delta will be adjusted dynamically.
*/ 
{
 public:
  static const int no_factor = -1;
  static const int dynamic_factor = -2;
  static const int special_factor = -3;
  int relevant_factor;
   // either the value is one out of { no_factor, dynamic_factor, special_factor }
   // or it is indexing a static factor (in the static factor base)

 private:
  CTinyFBsizetypeVector *Relation_sparse;
  myBitString *Relation_dense;
  mpz_t Delta; // Delta-value, that is the square-root part of the relation [a1*a2*a2 * (Delta)^2 = 1 mod n]
  std::istream *pDynamicRelations_from_file; // pointer to an istream, which provides us with correction factors (Korrekturfaktoren) for combine methods

 public:
  static void (*seek_emergency_handler)(istream &, const streampos);

  static void seek_emergency_default_handler(istream &, const streampos);

  class ProvideDynamicRelationsStream
   {
    public:
     typedef istream *Pistream; 
    private:
     Pistream &pistream;
     bool reuse;
    public:
     // constructor acquires an istream (usually for pDynamicRelations_from_file)
     ProvideDynamicRelationsStream(Pistream &pistream_); // acquire an istream for pistream_ (or reuse existing one)
     // destructor releases an istream (usually pDynamicRelations_from_file)
     ~ProvideDynamicRelationsStream(); // release istream (if it was acquired and not reused)
   }; 

 public:
  inline CRelation() // constructor
   : CProvideHelperVariables(), relevant_factor(no_factor),
     Relation_sparse(new CTinyFBsizetypeVector), Relation_dense(NULL),
     pDynamicRelations_from_file(NULL), pMulticombineData(NULL)
    {
      mpz_init_set_ui(Delta,1);
    }

  inline ~CRelation() // destructor
    {
#if 0
      if (Relation_sparse) { delete Relation_sparse; Relation_sparse=NULL; }
      if (Relation_dense) { delete Relation_dense; Relation_dense=NULL; }
      mpz_clear(Delta);
      relevant_factor=no_factor;
#else
      delete Relation_sparse;
      delete Relation_dense;
      mpz_clear(Delta);
#endif
    }
  
  explicit CRelation(const signed int SievePos, short int HitCount=0); // constructor
  /* initialize "Deltas" with the value of Delta,
     fill in "Relationen" with the prime factors (containing odd exponents)
     of the reduced square.
   */

 public:
  int largest_factor_in_Relation() const
    /* (index) of biggest factor in relation */
    {
      if (Relation_sparse)
        if (Relation_sparse->empty()) return no_factor;
        else 
          {
            return Relation_sparse->last();
          }
      else 
        {
          return Relation_dense->last(1);
          /*
            int h = Relation_dense->last(1);
            if (h>=0) return h;
            else return no_factor;
          */
        }
    }

  int second_largest_factor_in_Relation() const
    /* (Index) of second-largest factor in relation */
    {
      if (Relation_sparse)
        if (Relation_sparse->size()>1)
          return (*Relation_sparse)[Relation_sparse->size()-2];
        else
          return no_factor;
      else // Relation ist dense
        {
          signed int vorletzter_Faktor=Relation_dense->last(1);
          if (vorletzter_Faktor>=0)
           return Relation_dense->prev(vorletzter_Faktor,1);
          else
           return no_factor;
        }
    }

  inline unsigned int SizeOfRelation() const
    {
      if (Relation_sparse) return Relation_sparse->size();
      else return Relation_dense->count(1);
    }

  inline void optisize(void)
    {
      // optimize (minimize) memory consumption of the relation
      // 1. if appropriate, implicitly convert dense<->sparse,
      // 2. minimize memory consumption inside the representation

      const int size=SizeOfRelation();
      const int largest=largest_factor_in_Relation();
      const int quotient = size ? largest/size : 0;

      if (Relation_sparse)
       {
         if (quotient<16-2) convert_Relation_to_dense();
         else Relation_sparse->optisize();
       }
      else // (Relation_dense)
       {
         if (quotient>16+2) convert_Relation_to_sparse();
         else Relation_dense->optisize();
       }
    }
  
  void convert_Relation_to_dense();
  void convert_Relation_to_sparse();

  inline bool empty() const
    {
      return (relevant_factor==no_factor);
      // if (Relation_sparse) return Relation_sparse->empty();
      // else return Relation_dense->last(1)==-1; 
    }
        
  void combine(const CRelation& GL2);

  // routines for speeding up multiple combine steps
 public:

  // TExponentArrayElement:
  // This type is used for multi-combine routines to collect the exponents
  // according to each element of the static factorbase.  It should be an
  // unsigned type which is able to store the maximum count of multi-combine
  // calls between a multi-combine init and a multi-combine exit call. 
  // However, if the type is too small, intermediate "multi-combine-exit"
  // calls will be automatically generated to avoid overflows.
  // "unsigned short int" is a little bit faster than "unsigned int" for
  // mmx-aware machines (as it needs less memory and thereby also provides
  // better cache effectiveness).
  typedef unsigned short int TExponentArrayElement;

  struct SMulticombineData
   {
     TExponentArrayElement ExponentArray[StaticFactorbaseSettings::MaxSize+64] __attribute__ ((aligned (16)));
      // initial value must be {0} (please assure in the corresponding ".cc-file" that all elements are zero by default!)
      // +64 because we want to allow fast routines to avoid costly range checks...
      // 16bit alignment for fastest access in case that MMX/SSE is used to process this array
     unsigned long int multi_combine_Counter; // (initial value must be 0!)

#if 1 && defined(ASM_SSE2)
 #ifdef DEBUG
  #warning "CRelation::SMulticombineData: using own new/delete to overcome alignment problems"
 #endif
     // too bad: __attribute__ ((aligned (16))) does not work, if memory is allocated using "new".
     // Solution: a tiny new-operator that uses int posix_memalign(void **memptr, size_t alignment, size_t size);
     // this is only a tiny, ugly, basic replacement implementation to get things to work!!

     static void* operator new(size_t size)
      {
        void *p;
#ifdef CYGWIN_COMPAT
        p = memalign(16,size);
        if ( p==NULL || (reinterpret_cast<int>(p)&0x0f) )
         {
           cout << "NULL or not aligned p: " << p << endl;
           throw std::bad_alloc();
         }
#else
        if (posix_memalign(&p,16,size)) throw std::bad_alloc();
#endif
        return p;
      }
     static void operator delete(void *p, size_t /* size */)
      {
        ::free(p);
      }
#endif

     SMulticombineData() : multi_combine_Counter(0)
      {
#if defined(DEBUG) && (defined(ASM_SSE) || defined(ASM_SSE2))
        if (reinterpret_cast<unsigned int>(&ExponentArray[0]) &0x0fU)
         {
           MARK; cerr << "not 16-byte aligned!" << endl;
           cout << ExponentArray << endl;
           cout << this << endl;
         }
#endif
        for (int i=0; i<StaticFactorbaseSettings::Size(); i++) ExponentArray[i]=0;
      }
   };

 protected:
  SMulticombineData *pMulticombineData;
   // pointer to a MulticombineData structure, that will be used inside the
   // multicombine routines. Allocating/Deallocating the Structure is done outside of this class!

 public:
  void set_MulticombineData(SMulticombineData *Data) { pMulticombineData=Data; }
  void invalidate_MulticombineData() { pMulticombineData=NULL; }
  void dispose_MulticombineData() { delete pMulticombineData; invalidate_MulticombineData(); }
  void use_MulticombineData_from(const CRelation &GL) { pMulticombineData=GL.pMulticombineData; }
  void multi_combine_init();
  void multi_combine_main(const CRelation& GL2);
  void multi_combine_exit();

 protected:
  inline void adjust_multi_combine()
   {
     // increase multi_combine_Counter and protect against overflow
     if (++pMulticombineData->multi_combine_Counter>=std::numeric_limits<TExponentArrayElement>::max())
      {
#ifdef VERBOSE
        cout << "intermediate multi_combine_exit() issued" << endl;
#endif
        --pMulticombineData->multi_combine_Counter;
        multi_combine_exit(); // do intermediate combine_exit to avoid overflow
        ++pMulticombineData->multi_combine_Counter;
      }
   }

 public:
  bool ComputeQuadraticCongruence() const;
  bool is_valid() const; // check, whether the relation is valid (congruency must be valid)
  static bool is_valid(istream &in); // check, whether the relation given at the current file & fileposition is valid. (congruency must be valid)
  inline void SanityCheck() const // check, whether the relation is valid (congruency must be valid, otherwise: abort program)
   {
     if (!is_valid()) exit(1);
   }
  static void SanityCheck(istream &in) // check, whether the relation given at the current file & fileposition is valid. (congruency must be valid, otherwise: abort program)
   {
     if (!CRelation::is_valid(in)) exit(1);
   }
  static void SanityCheckRelationsFile(const std::string FileName); // check, whether the relations inside file are valid. (congruency must be valid, otherwise: abort program)

 private:
  inline void swap(CRelation &GL2)
    { 
      // swap content of this relation with the content of GL2
      // (using reference semantic)
      std::swap(Relation_sparse,GL2.Relation_sparse);
      std::swap(Relation_dense,GL2.Relation_dense);
      std::swap(relevant_factor,GL2.relevant_factor);
      mpz_swap(Delta,GL2.Delta);
    }
 public:
  inline bool operator< (const CRelation &GL2) const
    {   
      return relevant_factor < GL2.relevant_factor;
    }
  streampos save(ostream &out, const CmpqsFactor factor, const short int HitCount=0) const;
  inline streampos save(ostream &out, const int i=1, const short int HitCount=0) const
   {
     CmpqsFactor DLP;
     DLP=i;
     return save(out,DLP,HitCount);
   }
  CmpqsFactor combine(istream &in, const streampos pos);
  CmpqsFactor multi_combine_main(istream &in, const streampos pos);
  CmpqsFactor combine(istream &in);
  CmpqsFactor multi_combine_main(istream &in);
  friend ostream& operator<< (ostream& ostr, const CRelation& GL);
  friend class StaticRelations;
  friend class SpecialRelations;
  friend class Cprocess_clients;
};


ostream& operator<< (ostream& ostr, const CRelation& GL);



class TDynamicFactorRelation
{
  // If factor is not too big, then we can sieve with it as a dynamic factor.
  // At the moment dynamic factors are restricted to "int".
private:
  static void append_DynamicFactor_for_sieving(const int DynFac);
public:
  int factor; // factor is given as an integer value
  streampos fpos;

  TDynamicFactorRelation() : factor(0), fpos(-1) { }

  inline bool operator() (const TDynamicFactorRelation &t1, const TDynamicFactorRelation &t2) const
    { 
      // for STL: set
      return t1.factor<t2.factor;
    }
  inline int operator() (const TDynamicFactorRelation &t2) const
    { 
      // for GNU extension: hash_set
      return t2.factor;
    }
  inline bool operator== (const TDynamicFactorRelation &t2) const
    { 
      return factor==t2.factor;
    }
  inline bool sieveable() const
   {
     // are we allowed to sieve with this Large Prime?
     return factor<DynamicFactor_SievingThreshold;
   }
  inline void append_for_sieving() const
   {
#ifndef IS_SERVER
     // since the server itself does no sieving at all,
     // this is done only for clients and standalone versions:
     if (sieveable()) append_DynamicFactor_for_sieving(factor);
#endif
   }
};


struct TSieve_Delta
{
  int delta;
  int factor;
  inline bool operator< (const TSieve_Delta &v) const
    {
#if defined (USE_FIBHEAP)
      return (delta < v.delta);
#else
      return (delta > v.delta);
#endif
      // delta > v.delta , if queue is used
      // delta < v.delta , if fibonacci-heap is used
    }
};


#endif

---