CRelation.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <stdexcept>
#include <cmath>
#include <gmp.h>
#include "utils.H"

#include "qsieve.H"

using namespace std;

extern mpz_t n;


ostream& operator<< (ostream& ostr, const CRelation& GL)
{
  if (GL.Relation_dense) ostr << "dense"; else ostr << "sparse";
  ostr << "1=";
  
  char *str;
  str = new char [mpz_sizeinbase(GL.Delta,10)+2]; mpz_get_str(str, 10, GL.Delta);
  ostr << str << "^2";
  delete [] str; 
  
  if (GL.Relation_dense)
    { 
      for (int i=GL.Relation_dense->first(); i>=0; i=GL.Relation_dense->next(i))
        ostr << "*p" << i;
    }
  else
    {
      const CTinyFBsizetypeVector& R = *GL.Relation_sparse;
      for (int pos=0; pos<R.size(); ++pos)
        ostr << "*p" << R[pos];
    }
  
  return ostr;
}


void CRelation::convert_Relation_to_dense()
{
  //cout << "Converting relation from sparse to dense representation." << endl;
  if (Relation_dense)
    { cerr << "Relation is already dense!" << endl; return; }
  Relation_dense = new myBitString;
  for (int pos=Relation_sparse->size()-1; pos>=0; --pos)
    Relation_dense->set((*Relation_sparse)[pos]);
  //Relation_sparse->erase(Relation_sparse->begin(),Relation_sparse->end());
  delete Relation_sparse; Relation_sparse=NULL; // delete sparse representation
}

void CRelation::convert_Relation_to_sparse()
{
  //cout << "Converting relation from dense to sparse representation." << endl;
  if (Relation_sparse)
    { cerr << "Relation is already sparse!" << endl; return; }
  Relation_sparse = new CTinyFBsizetypeVector(Relation_dense->count(1));
  for (int i=Relation_dense->first(); i>=0; i=Relation_dense->next(i))
    Relation_sparse->fast_append(i);
  //Relation_dense->clear();
  delete Relation_dense; Relation_dense=NULL; // delete dense representation
}


void CRelation::combine(const CRelation& GL2)
  // GL2 as reference, because this more efficient
  // compute symmetric difference
{
  //cout << "combining relations..." << endl;
  //cout << "Gl.1: "; SanityCheck();
  //cout << "Gl.2: "; GL2.SanityCheck();
  
  /*
    Compute symmetric difference of the two relation sets.
    Multiply common part (intersection) into square root Delta.
  */
  mpz_set_ui(x,1); // collect factors in x until x>n (to speedup computation)

  if ( GL2.Relation_dense ) // GL2.Relation is dense
    {
      if (!Relation_dense) // if relation isn't dense already,
        convert_Relation_to_dense(); // make it dense...
      
      // PrimeNumbers[0] == -1; we can ignore the sign!
      //  -> loop until i>0 instead of i>=0
      //  -> we can treat all numbers as unsigned, thereby optimizing things...
      for (int i=GL2.Relation_dense->last(); i>0; i=GL2.Relation_dense->prev(i))
        if (Relation_dense->test(i))
         {
           mpz_mul_ui(x,x,PrimeNumbers[i]);
           if (mpz_cmp(x,n)>=0)
            {
              mpz_mod(x,x,n); mpz_mul(Delta,Delta,x);
              mpz_mod(Delta,Delta,n); mpz_set_ui(x,1);
            }
         }
      
      (*Relation_dense)._xor(*GL2.Relation_dense); // symmetric difference
    }
  else /* GL2.Relation is sparse */
    if (Relation_dense) /* and Relation is dense */
      { 
        // compute symmetric difference in dense representation
        // (whereas the other vector remains sparse)
        for (int pos=0; pos<GL2.Relation_sparse->size(); ++pos)
          {
            if (Relation_dense->test((*GL2.Relation_sparse)[pos]))
              if ((*GL2.Relation_sparse)[pos]>0)
                {
                  mpz_mul_ui(x,x,PrimeNumbers[(*GL2.Relation_sparse)[pos]]);
                  if (mpz_cmp(x,n)>=0)
                    {
                      mpz_mod(x,x,n); mpz_mul(Delta,Delta,x);
                      mpz_mod(Delta,Delta,n); mpz_set_ui(x,1);
                    }
                } //else mpz_neg(Delta,Delta);
            Relation_dense->invert((*GL2.Relation_sparse)[pos]);
          }
      }
    else
      {
        // compute symmetric difference in sparse representation
        int p1 = 0;
        int p2 = 0;
        const CTinyFBsizetypeVector& R = *Relation_sparse;
        const CTinyFBsizetypeVector& R2 = *GL2.Relation_sparse;
        CTinyFBsizetypeVector neueMenge(R.size()+R2.size());
        //neueMenge.reset();
        while (p2<R2.size() && p1<R.size())
          {
            while (R2[p2]<R[p1])
              {
                neueMenge.fast_append(R2[p2]); ++p2;
                if (p2==R2.size()) goto done;
              }
            while (R[p1]<R2[p2])
              {
                neueMenge.fast_append(R[p1]); ++p1;
                if (p1==R.size()) goto done;
              }
            while (R2[p2]==R[p1])
              {
                if (R[p1]>0)
                  {
                    mpz_mul_ui(x,x,PrimeNumbers[R[p1]]);
                    if (mpz_cmp(x,n)>=0)
                      {
                        mpz_mod(x,x,n); mpz_mul(Delta,Delta,x);
                        mpz_mod(Delta,Delta,n); mpz_set_ui(x,1);
                      }
                  } //else mpz_neg(Delta,Delta);
                ++p1; ++p2;
                if (p2==R2.size()) goto done;
                if (p1==R.size()) goto done;      
              }
          }
      done:
        while (p2!=R2.size()) { neueMenge.fast_append(R2[p2]); ++p2; }
        while (p1!=R.size()) { neueMenge.fast_append(R[p1]); ++p1; }
        Relation_sparse->copy_from(neueMenge); // take over new set
        //optisize(); // optimize memory usage (sparse<->dense, etc.)
      }
  
  // x contains (possibly) a remaining value, which has to be multiplied into Delta:
  mpz_mul(Delta,Delta,x); mpz_mod(Delta,Delta,n);
  
  // don't forget: the new delta has to contain the delta of the other
  // relation as well!
  mpz_mul(Delta,Delta,GL2.Delta); mpz_mod(Delta,Delta,n);
  
  relevant_factor=largest_factor_in_Relation();
  //cout << "combine finished." << endl;
}


// The following set of methods greatly improve multiple combination of relations.
//  idea: the update/combine of the exponents of the factors can
//   deferred up to the time, where the result is needed.
//   (In most situations there is really no need to this simultaneously.)

void CRelation::multi_combine_init()
{
  // initialize vector of exponents.
  // hint: If you can assure, that the last call of multi_combine_...
  //       was in fact a multi_combine_exit(), then you
  //       may skip calling multi_combine_init().
  if (pMulticombineData==NULL)
   {
     MARK; cerr << "Error: pMulticombineData invalid!" << endl;
     exit(1);
   }
  if (!pMulticombineData->multi_combine_Counter) return;
  pMulticombineData->multi_combine_Counter=0;
  for (int i=0; i<StaticFactorbaseSettings::Size(); i++) pMulticombineData->ExponentArray[i]=0;
}

void CRelation::multi_combine_exit()
{
  if (!pMulticombineData->multi_combine_Counter) return; // nothing to do...
#ifdef VERBOSE
  cout << "Multicombine_exit after " << pMulticombineData->multi_combine_Counter
       << " calls." << endl;
#endif

  if (pMulticombineData->multi_combine_Counter>std::numeric_limits<TExponentArrayElement>::max())
   {
     MARK;
     cerr << "Elements of CRelation::ExponentArray may have exceeded limits," << endl
          << "you have to increase the type of CRelation::TExponentArrayElement" << endl
          << "and recompile the program to assure correct computation." << endl
          << "(Data computed so far are still correct and can be reused afterwards.)" << endl;
     exit(8);
   }

  optisize(); // optimize Size of vector (sparse<->dense, etc.)

  // multiply factors weighted by their exponents into Delta:
  mpz_set_ui(y,1); // collect product in y
  pMulticombineData->ExponentArray[0]=0; // leave "prime number" "-1" out of consideration

  /* 
     implemented idea:
     sort vector of exponents descending by exponents,
     evaluate the product analogously to the Horner scheme.
     -> sorting method: bucketsort. ("Flag[i]" defines, whether an exponent is associated to i.)
     All bases belonging to an exponent e will be multiplied together in "Basis[e]".
  */

#define MY_BASIS_WORKAROUND
#if !defined(__GNUC__) || !defined(MY_BASIS_WORKAROUND)
  // I got strange errors in server mode when using this:
  // glibc 2.5 raises an assertion fault...
  // probably a racing condition between threads.
  mpz_t* const Basis = new mpz_t[StaticFactorbaseSettings::Size()];
  bool* const Flag = new bool[StaticFactorbaseSettings::Size()];
#else
  #warning "enabling workaround for possible Basis assertion failure"
  // this code is less standard conformant since it uses variable-length arrays,
  // but it seems to circumvent the above problem and may be slightly faster...
  __extension__ mpz_t Basis[StaticFactorbaseSettings::Size()];
  __extension__ bool Flag[StaticFactorbaseSettings::Size()];
#endif  
  for (int i=0; i<StaticFactorbaseSettings::Size(); ++i) Flag[i]=false; // initialize

  for (int i=1; i<StaticFactorbaseSettings::Size(); ++i)
   if (pMulticombineData->ExponentArray[i])
    {
      const unsigned long int e=pMulticombineData->ExponentArray[i];
      pMulticombineData->ExponentArray[i]=0; // reset exponent

      if (e>=StaticFactorbaseSettings::Size())
       {
         cout << "Unexpected exponent " << e << "! Use normal exponentiation..." << endl;
         mpz_set_ui(y,PrimeNumbers[i]); mpz_powm_ui(y,y,e,n);
         mpz_mul(Delta,Delta,y);
         continue;
       }

      if (Flag[e])
       {
         // expand Basis
         mpz_mul_ui(Basis[e],Basis[e],PrimeNumbers[i]);
         if (mpz_cmp(Basis[e],n)>=0)
          { 
            // reduce mod n
            mpz_mod(Basis[e],Basis[e],n);
          }
       }
      else 
       { 
         // initialize Basis
         Flag[e]=true;
         mpz_init_set_ui(Basis[e],PrimeNumbers[i]);
       }
    }

 /*
  Evaluate product analogously to Horner scheme,
  example: a^7 * b^4 * c^3 * d 
           = a^(7-4) * (a*b)^(4-3) * (a*b*c)^(3-1) * (a*b*c*d)
 */     

  unsigned int e2 = StaticFactorbaseSettings::Size()-1;
  Flag[0]=true; // be cautious for next loop...
  while (!Flag[e2]) --e2; // get highest set exponent

  mpz_set_ui(y,1);
  if (e2>1)
   for (unsigned int e=e2-1; e>0; --e)
    if (Flag[e])
     {
       mpz_mul(Basis[e],Basis[e],Basis[e2]);
       if (mpz_cmp(Basis[e],n)>=0)
        {
          // reduce mod n
          mpz_mod(Basis[e],Basis[e],n);
        }
       if (e2-e>1) mpz_powm_ui(Basis[e2],Basis[e2],e2-e,n);
       mpz_mul(y,y,Basis[e2]);
       if (mpz_cmp(y,n)>=0)
        {
          // reduce mod n
          mpz_mod(y,y,n);
        }
       mpz_clear(Basis[e2]); // do not forget: free this mpz-number
       e2=e;
     }

  if (e2>1) mpz_powm_ui(Basis[e2],Basis[e2],e2,n);
  if (e2>0)
   {
      mpz_mul(y,y,Basis[e2]); mpz_mod(y,y,n);
      mpz_clear(Basis[e2]); // do not forget: free this mpz number
   }

#if !defined(__GNUC__) || !defined(MY_BASIS_WORKAROUND)
  delete [] Basis; delete [] Flag;
#endif
  
  mpz_mul(Delta,Delta,y); mpz_mod(Delta,Delta,n);
  pMulticombineData->multi_combine_Counter=0;
#ifdef VERBOSE
  cout << "Multi_combine_exit finished." << endl;
#endif
}

void CRelation::multi_combine_main(const CRelation& GL2)
  // get GL2 via reference (to be more efficient)
  // compute symmetric difference
{
  //cout << "multi-combine relations..." << endl;
  adjust_multi_combine(); // keep aware of number of calls

  /*
    Compute symmetric difference of the two relation sets.
    Memorize common part (intersection) by adding it to ExponentArray.
  */

  if ( GL2.Relation_dense ) /* GL2.Relation is dense */
    {
      if (!Relation_dense) // relation is not dense
        convert_Relation_to_dense(); // make it dense!

      /* Zunächst werden die Bits der beiden Relationen
         (temporär!) komponentenweise durch "and" verknüpft.
         Dadurch werden die Primzahlen ermittelt, die verknüpft Quadrate
         ergeben. Diese können dann im Exponentenfeld berücksichtigt werden.
         In gewisser Weise wird die Bitstrings komponentenweise wie
         mit einem "Halbaddierer" verknüpft.
          Additionsergebnis (durch xor) überschreibt anschließend den
          Bitstring, während der Übertrag im Exponentenfeld addiert wird.
       */
      Relation_dense->test_and_add_carry(*GL2.Relation_dense,pMulticombineData->ExponentArray);
      (*Relation_dense)._xor(*GL2.Relation_dense); // compute symmetric difference
    }
  else // GL2.Relation is sparse
    if (Relation_dense) // and Relation is dense
      { // compute symmetric difference in dense representation,
        // whereas the other vector is sparse
        const CTinyFBsizetypeVector& R2 = *GL2.Relation_sparse;
        for (int pos=0; pos<R2.size(); ++pos)
          {
            if (Relation_dense->test_and_invert(R2[pos])) ++pMulticombineData->ExponentArray[R2[pos]];
          }
      }
    else
      { // compute symmetric difference in sparse representation
        int p1 = 0;
        int p2 = 0;
        const CTinyFBsizetypeVector& R  = *Relation_sparse;
        const CTinyFBsizetypeVector& R2 = *GL2.Relation_sparse;
        CTinyFBsizetypeVector neueMenge(R.size()+R2.size());
        //neueMenge.reset();
        while (p2<R2.size() && p1<R.size())
          {
            while (R2[p2]<R[p1])
              {
                neueMenge.fast_append(R2[p2]); ++p2;
                if (p2==R2.size()) goto done;
              }
            while (R[p1]<R2[p2])
              {
                neueMenge.fast_append(R[p1]); ++p1;
                if (p1==R.size()) goto done;
              }
            while (R2[p2]==R[p1])
              {
                ++pMulticombineData->ExponentArray[R[p1]];
                ++p1; ++p2;
                if (p2==R2.size()) goto done;
                if (p1==R.size()) goto done;      
              }
          }
      done:
        while (p2!=R2.size()) { neueMenge.fast_append(R2[p2]); ++p2; }
        while (p1!=R.size()) { neueMenge.fast_append(R[p1]); ++p1; }
        Relation_sparse->copy_from(neueMenge); // take over new set
      }

  // don't forget: the new delta has to contain the delta of the other
  // relation as well!
  mpz_mul(Delta,Delta,GL2.Delta); mpz_mod(Delta,Delta,n);
  
  //-------------------------------------------------------------------
  // IMPORTANT: After the last call of multi_combine_main it is necessary
  // to call a closing multi_combine_exit!
  // Otherwise the values in "Exponentenfeld" would be left out!
  //-------------------------------------------------------------------

  relevant_factor=largest_factor_in_Relation();
  //cout << "finished combine." << endl;
}


// default action, when an invalid seek occurs in one of the combine methods
void CRelation::seek_emergency_default_handler(istream &, const streampos)
{
  throw ios_base::failure("invalid seek position on istream");
}

// a function pointer, that gets called whenever combine tries to seek to an invalid position
void (*CRelation::seek_emergency_handler)(istream &, const streampos) = CRelation::seek_emergency_default_handler;

CmpqsFactor CRelation::combine(istream &in, const streampos pos)
{
  const streampos merkfpos=in.tellg(); // memorize fileposition
  in.seekg(pos); if (in.peek()==EOF) seek_emergency_handler(in,pos);
  CmpqsFactor f=combine(in);
  in.seekg(merkfpos); // restore memorized fileposition
  return f;
}

CmpqsFactor CRelation::multi_combine_main(istream &in, const streampos pos)
{
  const streampos merkfpos=in.tellg(); // memorize fileposition
  in.seekg(pos); if (in.peek()==EOF) seek_emergency_handler(in,pos);
  CmpqsFactor f=multi_combine_main(in);
  in.seekg(merkfpos); // restore memorized fileposition
  return f;
}


int CStreamDecoder::GetValue()
{
  if (ahead_value>0) { int w=ahead_value; ahead_value=0; return w; }

  char c;
  in >> c;
  if ( c>='A' && c<='Z' ) return c-('A'-1); // A=1,...,Z=26
  if ( c>='g' && c<='z' )
   {
     int h = (c-'g');
     ahead_value=(h%5)+1; // next "ahead read" value
     return (h/5)+1; 
   }
  if (c=='0') return 0; // as we have no leading zeros, '0' means 0.
  if (c=='$') return -1;
  if (c=='%') return -2;

  if (c==' ')
   {
     MARK;
     cerr << "GetValue: Read a space? How to handle this?" << endl;
     in.setstate(std::ios::failbit);
     throw std::runtime_error("GetValue: read invalid character");
   }

  // otherwise the memorized number:
  int w;
  in.unget();
  in >> w;
  return w;
}

void CStreamEncoder::PutValue(const int w)
{
  if (w<=5 && w>=1 && prev_value>0) // then output coded value
    {
      out << char( (prev_value-1)*5 + w-1 + 'g' );
      prev_value=0;
      return;
    }

  if (prev_value>0)
   {
     out << char(prev_value+('A'-1)); // output value
     prev_value=0;
   }

  if (w>26) 
   {
     out << w << " "; // output value without encoding
     return;
   }
  if (w>0)
   {
     if (w<=4) { prev_value=w; return; } // wait till next time...
     out << char(w+('A'-1)); // output value coded as char
     return;
   }
  switch (w)
   {
     case -2: out << "%"; break;
     case -1: out << "$"; break;
     case  0: out << "0"; break;
     default:
      MARK;
      cerr << "PutValue: don't know how to handle this..." << endl;
      cerr << "w=" << w << endl;
      out.setstate(std::ios::failbit);
      throw std::invalid_argument("PutValue: invalid value");
   }
}

bool CRelation::is_valid() const
{
  mpz_t x,y;
  mpz_init(x);
  mpz_init_set_ui(y,1);
  if (Relation_sparse)
    {
      for (int i=0; i<Relation_sparse->size(); ++i)
        {
          mpz_mul_si(y,y,PrimeNumbers[(*Relation_sparse)[i]]);
          mpz_mod(y,y,n);
        }
    }
  else
    {
      for (int i=Relation_dense->first(); i>=0; i=Relation_dense->next(i))
        {
          mpz_mul_si(y,y,PrimeNumbers[i]);
          mpz_mod(y,y,n);
        }
    }
  
  mpz_invert(x,Delta,n); mpz_powm_ui(x,x,2,n);

#ifdef VERBOSE
  cout << "checking for congruency: ";
  cout << x << " = " << y << " (mod N)" << endl;
#endif
  if (mpz_cmp(x,y)!=0) 
    {
#ifdef VERBOSE_WARN
      MARK; cout << "Relation is invalid: Values are NOT congruent!" << endl;
#endif
      mpz_clear(x); mpz_clear(y);
      return false;
    }
 mpz_clear(x); mpz_clear(y);   
 return true;
}


#include <fstream>
#include <set>
struct myintcompare_
{
  inline bool operator() (const int i1, const int i2) const
   {
     return i1 < i2;
   }
};
typedef std::set<int, myintcompare_> TMyIntSet;

bool CRelation::is_valid(istream &in)
{
  static TMyIntSet InductiveDynamicFactorSet; // set of already proven dynamic factors

  CStreamDecoder enc_in(in);
  bool retval = false;

  string s;
  CmpqsFactor MainFactor;

  mpz_t x,y,Delta;
  mpz_init_set_ui(x,1);
  mpz_init_set_ui(y,1);
  mpz_init_set_ui(Delta,1);
 
  // save file flags
  const ios::fmtflags oldFlags = in.flags(); // memorize old flags of the stream
  in.unsetf(ios::hex); in.setf(ios::dec); // enable decimal output (notice: because of "correction factors" there could be a recursive call!)
  in.unsetf(ios::showbase); // without any leading prefix
 
  in >> s; // starting symbol "G"
  if (s != "G")
    {
#ifdef VERBOSE_WARN
      MARK;
      cerr << "error: wrong position while reading relation" << endl;
#endif
      retval=false; goto done;
    }
  in >> MainFactor; // read factor (dynamic-factor or special-factor of relation or 1)
  in >> s; // delta
  mpz_set_str(Delta, s.c_str(), mpzbase_f);
  mpz_mod(Delta, Delta, n);

  // now we start to work in hexadecimal notation:
  in.unsetf(ios::dec); in.setf(ios::hex); // hexadecimal output
  in.unsetf(ios::showbase); // without "0x"-marker
  
  signed int distance; // hint: we work with differences instead of the values itself

  {
    const size_t sizeinbase_of_n = mpz_sizeinbase(n,2);
    mpz_set_ui(x,1); // for speedup: collect factors in x as long as x<n

    distance=enc_in.GetValue();
    int index=distance;
    while (distance >= 0)
     {
       if (index>=StaticFactorbaseSettings::Size())
        {
  #ifdef VERBOSE_INFO
          MARK;
  #endif
  #ifdef VERBOSE_WARN
        cerr << "static factor outside the factorbase (index " << index << ")" << endl;
  #endif
        retval=false; goto done;
        }

       mpz_mul_si(x,x,PrimeNumbers[index]);
       if (mpz_sizeinbase(x,2)>sizeinbase_of_n)
        {
          mpz_mod(x,x,n); mpz_mul(y,y,x);
          mpz_mod(y,y,n); mpz_set_ui(x,1);
        }
       distance=enc_in.GetValue(); index+=distance;
     }
    mpz_mul(y,y,x); mpz_mod(y,y,n);
  }
  
  if (distance == -2) // "correction factors" follow
    {
      /* Korrekturfaktoren sind dynamische Faktoren, die zum Zeitpunkt der Generierung
         der gelesenen Relation nicht aufgelöst wurden. Sie müssen daher jetzt verknüpft werden. */
      /* correction factors are dynamic factors, that weren't resolved at the time of generating
         this relation. They must be considered and combined now! */
      int Korrekturfaktor;
      while ( (Korrekturfaktor=enc_in.GetValue())>0 )
       {
          const TMyIntSet::const_iterator p = InductiveDynamicFactorSet.find(Korrekturfaktor);
          if (p==InductiveDynamicFactorSet.end())
            {
#ifdef VERBOSE_INFO
              MARK;
#endif
#ifdef VERBOSE_WARN
              cerr << "Unknown correction value (Korrekturfaktor): " << Korrekturfaktor << endl;
#endif
              retval=false; goto done;
            }
          mpz_mul_ui(y,y,Korrekturfaktor); mpz_mod(y,y,n);
        }
    }

  if (MainFactor.LP1()) mpz_mul_si(y,y,MainFactor.LP1());
  if (MainFactor.LP2()) mpz_mul_si(y,y,MainFactor.LP2());
  mpz_mod(y,y,n);

  mpz_invert(x,Delta,n); mpz_powm_ui(x,x,2,n);
  
#ifdef VERBOSE
  cout << "checking for congruency: ";
  cout << x << " = " << y << " (mod N)" << endl;
#endif
  if (mpz_cmp(x,y)!=0) 
    {
#ifdef VERBOSE_WARN
      MARK; cout << "Relation is invalid: Values are NOT congruent!" << endl;
#endif
      retval=false; goto done;
    }

  if (MainFactor.IsTypeOf(CmpqsFactor::single_large_prime))
   {
     // now, since this relation is correct, other relations that
     // rely on this dynamic factor (as Korrekturfaktor)
     // can safely use this factor without destroying their correctness
     InductiveDynamicFactorSet.insert(MainFactor.int_value());
   }

  retval=true; // relation is okay
done:
  mpz_clear(Delta); mpz_clear(y); mpz_clear(x);
  in.flags(oldFlags); // restore memorized stream flags
  return retval;
}

void CRelation::SanityCheckRelationsFile(const std::string FileName)
{
#ifdef VERBOSE_NOTICE
  cout << "Performing sanity check for " << FileName << endl;
#endif
  std::ifstream MyFile(FileName.c_str());
  if (!MyFile)
   {
     cerr << "problems to open \"" << FileName << "\"!" << endl;
     exit(1);
   }

  MyFile.seekg(0,ios::beg);
  unsigned int count = 0, line=0, invalid=0, obsolete=0;
  while (MyFile)
   {
     while (MyFile.peek()=='U') { ++line; ++obsolete; MyFile.ignore(1000000,'\n'); }
     while (MyFile.peek()=='G')
      {
        ++count; ++line;
        if (!is_valid(MyFile))
         {
           ++invalid;
           cerr << "invalid Relation in " << FileName << ", line " << line << endl;
           MyFile.ignore(1000000,'\n');
         } else MyFile.ignore(1,'\n');
       #ifdef VERBOSE_NOTICE
        if (count%500==0) cout << count << " relations validated.\r" << flush;
       #endif
      }
     if (MyFile.peek()!=EOF && MyFile.peek()!='U')
      {
        ++line; ++invalid;
        cerr << "Garbage in " << FileName << ", line " << line << endl;
        cerr << "Giving up!" << endl;
        break;
        //MyFile.ignore(1000000,'\n');
      }
     // it seems that gcc-3.2.3 has some problems detecting eof,
     // I hope this helps...
     MyFile.get();
     if (MyFile.eof()) break;
     MyFile.unget();
   }
#ifdef VERBOSE_NOTICE
  cout << line << " lines scanned.       " << endl;
  if (obsolete) cout << obsolete << " obsolete relations skipped." << endl;
  cout << count << " relations validated." << endl;
#endif
  if (invalid)
   {
     cout << invalid << " INVALID relations detected in " << FileName << endl;
     exit(1);
   }
}

---