fermat.cc

Go to the documentation of this file.

void fermat_like_method()
{
  /* Fermat-ähnliche Methode:
     Sei N=p*q (hier ausnahmsweise p,q nichtnotwendigerweise prim)
     und sei M=m1*m2 ein Multiplikator, so dass m1*q ungefähr m2*p entspricht.
     Dann ist (m1*q)*(m2*p)=(m1*m2)*(q*p)=M*N.

     Wir definieren R=floor(sqrt(M*N)).
     Der Ansatz ist nun, ein delta zu finden, so dass Q := (R+delta)^2 - M*N
     eine Quadratzahl ergibt, denn dann folgt (R+delta)^2 - Q = 0 (mod N),
     was gleichbedeutend ist zu (R+delta+sqrt(Q))*(R+delta-sqrt(Q)) = 0 (mod N).
     Mit gcd(R+delta+sqrt(Q),N) wäre dann ein Teiler von N gefunden...
       
     Diese Methode ist zwar nicht besonders effizient, da sie eher einer
     Probedivision entspricht; sie hilft aber, solche Zahlen besonders
     schnell zu faktorisieren, die auf "allzu naive" Weise als
     "besonders schwer" zu faktorisierende Zahlen konstruiert sind.
     So können zum Beispiel Produkte der Form
       N:= p*nextprime(10*p+"kleiner Offset")
     effizient zerlegt werden.
  */
#ifdef VERBOSE_NOTICE
  cout << "starting Fermat factorization method" << endl;
#endif
  if (mpz_sizeinbase(n,10)<50)
   {
     cout << "number too small, it's not worthwhile do proceed; skipping fermat..." << endl;
     return;
   }
  
  const unsigned int Witnesses = 4;
  const unsigned int Witness[Witnesses] = { 23*41*59, 35*29*37, 31*43*47, 61*67 };
  const unsigned int Maxbits = 65536;
  unsigned char squareflags[Maxbits]; 

  // initialize bitfield for witnessing non-squares via quadratic residues
  for (unsigned int j=0; j<Maxbits; ++j) squareflags[j]=0;
  for (unsigned int i=0; i<Witnesses; ++i)
   for (unsigned int j=0; j<Witness[i]; ++j) squareflags[(j*j)%Witness[i]] |= 1<<i;

  mpz_t c,r,h;
  mpz_init(c); mpz_init(r); mpz_init(h);

  const unsigned int k_max = 0x100; // or 0x040 (to be less exhaustive)
  for (unsigned int k=0; k<k_max; ++k)
    {
#ifdef REACT_ON_SIGUSR
      if (USRSignalHandler.got_SIGUSR1()) break;
      if (USRSignalHandler.got_SIGUSR2()) break;
#endif
      mpz_set_ui(c,2*2*3*5*7*11*13*17*19); mpz_pow_ui(c,c,4);
      unsigned int m = 1;
      if (k&0x001) m*=2;
      if (k&0x002) m*=3;
      if (k&0x004) m*=5;
      if (k&0x008) m*=7;
      if (k&0x010) m*=11;
      if (k&0x020) m*=13;
      if (k&0x040) m*=17;
      if (k&0x080) m*=19;
      mpz_mul_ui(c,c,m);

      // Idea behind these multipliers given by example:
      // Let t be a (small) number, then
      //
      //              t^0   t^1   t^2   t^3   t^4
      //  t^4 covers  --- , --- , --- , --- , --- 
      //              t^4   t^3   t^2   t^1   t^0
      //
      //  which is identical to the sequence
      //
      //    1     1     
      //   --- , --- , 1 , t^2 , t^4
      //   t^4   t^2
      //
      // because p*q*t^4 = p*(q*t^4) = (p*t)*(q*t^3) = (p*t^2)*(q*t^2) = ...
      //
      //                        1     1    1
      // Similarly t^5 covers  --- , --- , - , t , t^3, t^5
      //                       t^5   t^3   t

#if defined(VERBOSE_INFO)
      cout << "using Fermat multiplier " << c << "             \r" << flush;
#endif
      mpz_mul(c,c,n); mpz_sqrt(r,c); mpz_mul(h,r,r); mpz_sub(h,h,c);
      mpz_mul_ui(r,r,2);

      struct { unsigned int h,r; } ring[Witnesses];
      for (unsigned int i=0; i<Witnesses; ++i)
       {
         ring[i].h=mpz_remainder_ui(h,Witness[i]);
         ring[i].r=mpz_remainder_ui(r,Witness[i]);
         //cout << "Witness " << i << ": " << ring[i].h << ", " << ring[i].r << endl;
       }

      for (unsigned int delta=1; delta<50000; ++delta)
        {
          //if (delta%10000==0) cout << "delta=" << delta << "\r" << flush;
          mpz_add(h,h,r); mpz_add_ui(h,h,1); mpz_add_ui(r,r,2); // compute next square: h=(x+1)^2=x^2+2x+1=h+r+1 (using r=2x)
          for (unsigned int i=0; i<Witnesses; ++i)
           {
             ring[i].h+=ring[i].r+1; if (ring[i].h>=Witness[i]) ring[i].h-=Witness[i];
             ring[i].r+=2; if (ring[i].r>=Witness[i]) ring[i].r-=Witness[i];
           }

          if (!(squareflags[ring[0].h] & 0x01)) continue;
          if (!(squareflags[ring[1].h] & 0x02)) continue;
          if (!(squareflags[ring[2].h] & 0x04)) continue;
          if (!(squareflags[ring[3].h] & 0x08)) continue;
          //cout << "delta=" << delta << endl;          
          if (mpz_perfect_square_p(h))
            {
#if 1 || defined(VERBOSE_INFO)
              cout << endl << "delta=" << delta << ": Square detected..." << endl;
#endif
              mpz_div_ui(r,r,2); mpz_sqrt(h,h); mpz_add(h,h,r);
              mpz_gcd(h,h,n);
              if (mpz_cmp_ui(h,1)!=0) 
                {
                  mpz_divexact(n,n,h);
                  //cout << "t1=" << h << endl;
                  //cout << "t2=" << n << endl;
                  if (mpz_probab_prime_p(h,probab_prime_checks))
                    {
#if 0
 // only for debugging...
{
  mpz_t sum,N,x;
  mpz_init(sum); mpz_init(N); mpz_init(x);
  mpz_add(sum,n,h); mpz_div_ui(sum,sum,2); mpz_mul(N,n,h);
  mpz_sqrt(x,N); mpz_sub(sum,sum,x);
  cout << "phimat-value --> " << sum << endl;
  mpz_clear(x); mpz_clear(N); mpz_clear(sum); 
}
#endif

                      cout << h << " is factor." << endl;
                      std::ostringstream comment;
                      comment << " [fer]";
                      Factorization_to_file << MAL(h,comment) << flush;
                    }
                  else
                    {
                      fermat_like_method(); // try again...
                      mpz_swap(n,h);
                      if (mpz_probab_prime_p(h,probab_prime_checks))
                       {
                         cout << h << " is factor." << endl;
                         std::ostringstream comment;
                         comment << " [fer]";
                         Factorization_to_file << MAL(h,comment) << flush;
                       }
                      else
                       {
                         cout << h << " is a composite factor." << endl;
                         if (Potenztest(h)) Factorization_to_file << " [fer]" << flush;
                         else
                          {
                            std::ostringstream comment;
                            comment << " [fer] [composite]";
                            Factorization_to_file << MAL(h,comment) << flush;
                          }
                       }
                    }
                  if (!mpz_probab_prime_p(h,probab_prime_checks)) fermat_like_method(); // and try again once more...
                  k=k_max; break;
                }
            }
        }
    }
#ifdef VERBOSE_INFO
  cout << endl;
#endif
  mpz_clear(h); mpz_clear(r); mpz_clear(c);
}


#include "utils.H"
#include <vector>
#include <algorithm>

  class entry
  {
   private:
    mp_limb_t hashval;
    friend class phimahashvecs;
    static unsigned int iL;
    static mpz_t L;
   public:
    static mpz_t Base;
    entry(const unsigned int index)
     {
       if (index==iL+1)
        {
          // shortcut, inc case that items are inserted "in order"
          ++iL; mpz_mul(L,L,Base); mpz_mod(L,L,n);
        }
       else
        {
          mpz_powm_ui(L,Base,index,n);
        }
       hashval=mpz_getlimbn(L,0);
     }

    ~entry() { }
    const mpz_t& get_mpz(const unsigned int index) const
     {
       if (index==iL) return L;
       else
        {
          mpz_powm_ui(L,Base,index,n);
          return L;
        }
     }
    mp_limb_t get_hash() const { return hashval; }
    bool operator< (const entry &rhs) const
     {
       // important: we do only compare the hash values, not the multiple
       // precision numbers!
       return (hashval<rhs.hashval);
     }
    bool operator== (const entry &rhs) const
     {
       return (hashval==rhs.hashval);
     }
    signed int compare(const mp_limb_t rhs) const
     {
       // compare based on hash value
       return (hashval<rhs) ? -1 : (hashval>rhs) ? 1 : 0;
     }

  };

 class phimavec : private std::vector<entry>
 {
  private:
   friend class phimahashvecs;
   double d_linear_approx_multiplicator; 
  public:
   phimavec() : d_linear_approx_multiplicator(0.0) { }
   ~phimavec() { }
   void sort()
    {
      //cout << "sorting " << size() << " items." << endl;
      std::sort(begin(),end());
      //size_t s=size();
      erase(std::unique(begin(),end()),end()); // remove duplicate hash values
      //if (size()!=s) cout << s-size() << " duplicates removed." << endl;
      d_linear_approx_multiplicator = static_cast<double>(size())/(at(size()-1).get_hash()-at(0).get_hash());
    }
   bool found(const mpz_t &val) const
    {
#ifdef VERBOSE_INFO
      static unsigned int cs = 0;
#endif

      // binary search
      signed int L=0, R=size()-1;
      const mp_limb_t hash = mpz_getlimbn(val,0);

#if 1
      // linear approximation of the index
      // (uniform distribution -> much better than binary search!)
      const signed int d = static_cast<signed int>(static_cast<double>(hash) * d_linear_approx_multiplicator);
      //if ( d<L || d>R ) return false; // should be outside the interval (but there are rounding errors...)
      const signed int dd = 127;
      __builtin_prefetch(&operator[](d),0,1);
      if (d-dd>L && at(d-dd).get_hash()<hash) L=d-dd;
      if (d+dd<R && at(d+dd).get_hash()>hash) R=d+dd;
      //cout << "-->L: " << L << " R: " << R << " d: " << d << endl;
#endif

      while (R>=L)
       {
         //cout << "searching " << L << " - " << R << endl;
         const signed int pos = (L+R)>>1;
         __builtin_prefetch(&operator[]((L+pos-1)>>1),0,0);
         __builtin_prefetch(&operator[]((R+pos+1)>>1),0,0);
         const signed int vgl = operator[](pos).compare(hash);
#ifdef VERBOSE_INFO
         ++cs;
#endif
         if (vgl==0)
          {
#ifdef VERBOSE_INFO
            cout << cs << " comparisons performed so far" << endl; 
#endif
            return true;
          }
         if (vgl<0) L=pos+1; else R=pos-1;
       }
      //cout << "L: " << L << " R: " << R << " d: " << d << " size:" << size() << endl;

      return false; 
    }
 };


 class phimahashvecs
 {
  private:
   static const mp_size_t hashlimb = 1; // for mpz_getlimbn(<number>,hashlimb)
   static const unsigned int lookup_bins = 0x1000;
   static const mp_limb_t limbmask = lookup_bins-1;
   phimavec lookup[lookup_bins];
   size_t myIntervalsize;
   const static unsigned int redundancy = 8;
  public:
   phimahashvecs() : myIntervalsize(0)
    {
      entry::iL=0; mpz_init_set_ui(entry::L,1);
    }
   ~phimahashvecs() { }

   size_t size() const
    {
      size_t s = 0;
      for (unsigned int i=0; i<lookup_bins; ++i) s+=lookup[i].size();
      return s;
    }
   size_t max_size() const { return lookup[0].max_size(); }
 
   void insert(const unsigned int i)
    {
      const entry e(i);
      const mp_limb_t h=mpz_getlimbn(e.get_mpz(i),hashlimb) & limbmask;
      lookup[h].push_back(e);
    }
   void prepare(const size_t tablesize)
    {
      // for optimal RAM requirements we precalculate the number of items
      // in each lookup bin.

      cout << "optimizing memory layout" << endl;
      entry::iL=0; mpz_set_ui(entry::L,1);
      size_t bin_size[lookup_bins] = { 0 };
      for (unsigned int i=0; i<tablesize; ++i)
       {
         const mp_limb_t h=mpz_getlimbn(entry::L,hashlimb) & limbmask;
         ++bin_size[h];
         mpz_mul(entry::L,entry::L,entry::Base); mpz_mod(entry::L,entry::L,n);
       }
      entry::iL=0; mpz_set_ui(entry::L,1);
      for (unsigned int i=0; i<lookup_bins; ++i)
       {
         //cout << "binsize[" << i << "]=" << bin_size[i] << endl;
         lookup[i].reserve(bin_size[i]);
       }
 
      // and now actually insert the items
      cout << "inserting items" << endl;
      for (size_t i=0; i<tablesize; ++i) insert(i); // implicit constructor call

      myIntervalsize = size()-redundancy; // important! call it before duplicate hash values are erased!

#ifdef VERBOSE_INFO
      size_t s = 0;
      for (unsigned int i=0; i<lookup_bins; ++i)
       {
         lookup[i].sort();
         s+=lookup[i].size();
         cout << s << " items prepared.\r" << flush;
       }
      cout << endl;
#else
      for (unsigned int i=0; i<lookup_bins; ++i) lookup[i].sort();
#endif

    }

   size_t Intervalsize() const { return myIntervalsize; }

   bool found(const mpz_t &val) const
    {
      // first hash into lookup, then do binary search inside this bin
      const mp_limb_t hash = mpz_getlimbn(val,hashlimb) & limbmask;
      if (!lookup[hash].found(val)) return false;

      // now it gets more complicated:
      // a possible candidate has been found based on equal hash values;
      // perform up to <redundancy> additional checks to reduce the
      // risk of false positives
      //cout << "candidate found; checking..." << endl;
      bool retval=true;
      mpz_t x;
      mpz_init_set(x,val);
      for (unsigned int i=0; i<redundancy; ++i)
       {
         mpz_mul(x,x,entry::Base); mpz_mod(x,x,n);
         const mp_limb_t hash = mpz_getlimbn(x,hashlimb) & limbmask;
         retval=lookup[hash].found(x);
         //cout << "candidate " << ((retval) ? "passed" : "failed") << " test " << i << endl;
         if (!retval) break;
       }
     mpz_clear(x);
     return retval;
    }

   bool found(const mpz_t &val, unsigned int &the_index) const
    {
      // first hash into lookup, then do binary search inside this bin
      if (!found(val)) return false;

      // a candidate has been found with high probability
      // to determine "the_index", we perform a binary search!

      mpz_t x;
      mpz_init(x);
      size_t d = Intervalsize()-1;
      { 
        register unsigned int i = 0;
        while (d>>=1) ++i;
        d=1<<i;
      }
      cout << "candidate found; performing binary search for the_index" << endl;
      size_t hhh=0;
      while (d)
       {
         mpz_powm_ui(x,entry::Base,hhh+d,n);
         mpz_invert(x,x,n); mpz_mul(x,x,val); mpz_mod(x,x,n);
         if (found(x)) hhh+=d;
         d>>=1;
       }
      mpz_powm_ui(x,entry::Base,hhh,n);
      if (mpz_cmp(x,val)!=0)
       {
         cout << "binary search failed? try slowsearch..." << endl;
         mpz_t y; mpz_init(y); mpz_invert(y,entry::Base,n);
         while (hhh>0)
          {
            mpz_mul(x,x,y); mpz_mod(x,x,n); --hhh;
            if (mpz_cmp(x,val)==0) break;
          }
         mpz_clear(y);
       }
      the_index=hhh;
      cout << "the_index = " << the_index << endl;
      mpz_clear(x);
      return true;
    }
 };


unsigned int entry::iL = 0;
mpz_t entry::L, entry::Base;

extern unsigned long int allowed_memory_usage_KB(void);


class lambda_delta
{
 private:
   mpf_t SQRTN, DELTA;
   mpf_t fres, f1res;
   static const unsigned int bits = 2048;
 public:

  static void Delta_by_ratio(mpz_t &Delta, const mpz_t n, const double ratio_p=1.0, const double ratio_q=1.0)
   {
     mpf_set_default_prec(bits);
     mpf_t x,y;
     mpf_init(x); mpf_init(y);
     mpf_set_d(x,ratio_p); mpf_set_d(y,ratio_q);
     mpf_div(x,x,y); // full precision division (important!)
     mpf_sqrt(x,x);
     mpf_ui_div(y,1,x); mpf_add(x,x,y); mpf_sub_ui(x,x,2);
     mpf_set_z(y,n); mpf_sqrt(y,y); mpf_mul(x,x,y);
     mpz_set_f(Delta,x);
   }


  lambda_delta(const mpz_t n, const mpz_t delta)
   {
     mpf_set_default_prec(bits);
     mpf_init(SQRTN); mpf_set_z(SQRTN,n); mpf_sqrt(SQRTN,SQRTN); // mpf_floor(SQRTN,SQRTN);
     mpf_init(DELTA); mpf_set_z(DELTA,delta);
     mpf_init(fres); mpf_init(f1res);
   }
  ~lambda_delta()
   {
     mpf_clear(f1res); mpf_clear(fres);
     mpf_clear(DELTA); mpf_clear(SQRTN);
   }

  const mpf_t& f(const mpf_t &x)
   {
     mpf_t h;
     mpf_init(h);
     mpf_sqrt(fres,x); mpf_ui_div(h,1,fres); mpf_add(fres,fres,h);
     mpf_sub_ui(fres,fres,2);
     mpf_mul(fres,fres,SQRTN);
     mpf_sub(fres,fres,DELTA);
     mpf_clear(h);
     return fres;
   }
  const mpf_t& f1(const mpf_t &x)
   { 
     mpf_t h;
     mpf_init(h);
     mpf_sqrt(f1res,x); mpf_ui_div(f1res,1,f1res);
     mpf_pow_ui(h,f1res,3); mpf_sub(f1res,f1res,h);
     mpf_div_ui(f1res,f1res,2);
     mpf_mul(f1res,f1res,SQRTN);
     mpf_clear(h);
     return f1res;
   }
#if 0   
  void newton()
   {
     mpf_t x0,x1;
     mpf_init(x0);
     mpf_init_set_d(x1,1.1);
     do
      {
        mpf_set(x0,x1);
        mpf_div(x1,f(x0),f1(x0)); mpf_sub(x1,x0,x1);
        
      } while (mpf_eq(x1,x0,bits-64)==0);
     cout << "newton: "; mpf_out_str(NULL,10,50,x1); cout << endl;
     //mpf_sqrt(x0,x1); mpf_div(x0,SQRTN,x0);
     //mpz_t p_bound; mpz_init(p_bound); mpz_set_f(p_bound,x0);
     //cout << "N=p*q, p<" << p_bound << endl;
     //int i; std::cin >> i;
     mpf_clear(x1); mpf_clear(x0);
   }
#endif     
};


void phimat2()
{
  cout << "Starting experimental Phimat factoring method." << endl;
  cout << "Warning: This method runs very long except for special constructed numbers" << endl;
  cout << "of the form N=p*q; p,q prime and q=p+x, x small compared to p." << endl;

  mpz_init_set_ui(entry::Base,2);

  mpz_t V,x;
  mpz_init(x); mpz_sqrt(x,n); mpz_mul_ui(x,x,2);

  unsigned int step=2;
  mpz_t Delta,initial_Delta,stop_Delta;
  mpz_init_set_ui(Delta,0);
  mpz_init(stop_Delta);

  // if this small block is omitted, then default ratio 1:1 is used.
  // Otherwise calculate Delta based on the assumption, that
  // n=p*q and p_ratio*p = q_ratio*q  (approximately)
  lambda_delta::Delta_by_ratio(Delta,n,1.0,1.0);
  lambda_delta::Delta_by_ratio(stop_Delta,n,2.0,1.0);
  
  mpz_sub_ui(Delta,Delta,mpz_remainder_ui(Delta,3*32));
  mpz_init_set(initial_Delta,Delta);

  {
    unsigned int r = mpz_remainder_ui(n,32);
    unsigned int disp=0; // displacement
    if (r%4==1) { disp=2; step=4; }
    else if (r%8==3) { disp=4; step=8; }
    else if (r%8==7) { disp=0; step=8; }
    else { cerr << "unknown state " << r << endl; exit(1); }
    if (mpz_remainder_ui(n,6)==5) { step*=3; disp*=3; }
    cout << "phimat multiplier is " << step << endl;
    step/=2;

    r=mpz_remainder_ui(x,step);
    mpz_sub_ui(Delta,Delta,r); mpz_add_ui(Delta,Delta,disp);
    mpz_sub_ui(x,x,r); mpz_add_ui(x,x,disp);
  }

  mpz_init_set(V,n); mpz_add_ui(V,V,1); mpz_sub(V,V,x);
  mpz_powm(V,entry::Base,V,n);
  mpz_powm_ui(entry::Base,entry::Base,step,n);

  // remark:
  // the "Delta" value for phimat should be exactly twice the value
  // of the "Fermat"-Delta value for the same number.
  // This is because "phimat" tries to determine p+q, whereas
  // the "Fermat" procedure calculates implicitly (p+q)/2.
  // Both procedures terminate, when "Delta" has reached this value.
  // (Well, maybe phimat terminates earlier with a very low probability,
  // because the "congruency" 2^(N+1) == 2^(p+q) (mod N) is not necessarily
  // a (perfect) "identity" (and we do not know for certain, that p,q are the
  // only prime factors of N).


  phimahashvecs entries; // structure of indices (sorted by powers)

  {
    //const unsigned int tablesize = 0x400000;
    const unsigned int t1 = entries.max_size();
    unsigned int t2 = static_cast<size_t>((allowed_memory_usage_KB()/sizeof(entry)));
    if (t2>65536) t2-=32768;
    t2*=1024;
    const unsigned int tablesize = MIN(t1,t2);
    cout << "Creating table of " << tablesize << " items." << endl;
    cout << "This may take a while..." << endl;
    entries.prepare(tablesize);
    cout << "Okay, items are prepared. Now searching for factors..." << endl;
  }


  {
    // correction for initial Delta (in case it's not zero)
    mpz_t x;
    mpz_init(x);
    mpz_div_ui(x,initial_Delta,step);
    mpz_powm(x,entry::Base,x,n);
    mpz_invert(x,x,n); // we walk downwards...
    mpz_mod(x,x,n);
    mpz_mul(V,V,x); mpz_mod(V,V,n);
    mpz_clear(x);
  }

  mpz_t PowIntervalStep;
  mpz_init_set(PowIntervalStep,entry::Base);
  mpz_powm_ui(PowIntervalStep,PowIntervalStep,entries.Intervalsize(),n);
  mpz_invert(PowIntervalStep,PowIntervalStep,n); // we walk downwards...
  mpz_mod(PowIntervalStep,PowIntervalStep,n);

  mpz_t TableStep;
  mpz_init_set_ui(TableStep,entries.Intervalsize()); mpz_mul_ui(TableStep,TableStep,step);

  for (unsigned int i=0; /* endlessly */ ; ++i)
   {
     unsigned int index;
     if (entries.found(V,index))
      {
        mpz_t x,y,z,myDelta;
        mpz_init(x); mpz_init(y); mpz_init(z); mpz_init_set(myDelta,Delta);
        mpz_set_ui(x,index); mpz_mul_ui(x,x,step); mpz_add(myDelta,myDelta,x);
        cout << "Found a candidate at Delta=" << myDelta << endl;
        mpz_sqrt(x,n); mpz_mul_ui(x,x,2); mpz_add(x,x,myDelta);
        mpz_mul(y,x,x); mpz_mul_ui(z,n,4); mpz_sub(y,y,z);
        mpz_sqrt(y,y); mpz_add(y,y,x);
        mpz_gcd(z,y,n); mpz_mod(z,z,n);
        if (mpz_cmp_ui(z,1)>0)
         {
           cout << "looks good :-)" << endl;
           mpz_set(Delta,myDelta); // accepted... :-)

           while (mpz_cmp_ui(n,1)!=0)
            {
             if (mpz_probab_prime_p(z,probab_prime_checks))
              {
                cout << z << " is factor." << endl;
                std::ostringstream comment;
                comment << " [phimat]";
                Factorization_to_file << MAL(z,comment) << flush;
              }
             else
              {
                cout << z << " is a composite factor." << endl;
                if (Potenztest(z)) Factorization_to_file << " [phimat]" << flush;
                else
                 {
                   std::ostringstream comment;
                   comment << " [phimat] [composite]";
                   Factorization_to_file << MAL(z,comment) << flush;
                 }
              }
             mpz_divexact(n,n,z);
             mpz_set(z,n);
            }
         }
        mpz_clear(myDelta); mpz_clear(z); mpz_clear(y); mpz_clear(x);
        if (mpz_cmp_ui(n,1)==0) break;
      }
     mpz_add(Delta,Delta,TableStep);
     mpz_mul(V,V,PowIntervalStep); mpz_mod(V,V,n);
     if ((i&0x3ffff)==0)
      {
        if ((i&0x1ffffff)==0) // status report
         {
           mpz_t x,y,z;
           mpz_init(x); mpz_init(y); mpz_init(z);
           mpz_sqrt(x,n); mpz_mul_ui(x,x,2); mpz_add(x,x,Delta);
           mpz_mul(y,x,x); mpz_mul_ui(z,n,4); mpz_sub(y,y,z);
           mpz_sqrt(y,y); mpz_add(y,y,x);
           mpz_div_ui(y,y,2);
           mpz_div(x,n,y);
           cout << "new bound p<" << x << endl;
           mpz_clear(z); mpz_clear(y); mpz_clear(x);
           //cout << "stop_Delta: " << stop_Delta << endl;
           if (mpz_cmp(Delta,stop_Delta)>0) break;
         }
        cout << "Delta=" << Delta << '\r' << flush;
      }
   }

  // for validation reasons give the Delta for Fermat,
  // where the square could be found...
  mpz_div_ui(Delta,Delta,2);
  cout << "Fermat Delta would be " << Delta << endl;

  mpz_clear(TableStep); mpz_clear(PowIntervalStep);
  mpz_clear(initial_Delta); mpz_clear(stop_Delta); mpz_clear(Delta);
  mpz_clear(x); mpz_clear(V);

  mpz_clear(entry::Base);
}



void phimat(mpz_t n)
{
  mpz_t x,y;
  mpz_init(x); mpz_init(y);
  
  mpz_mul_ui(x,n,4); mpz_sqrt(x,x); // x=trunc(2*sqrt(n));
  mpz_add_ui(y,n,1); mpz_sub(y,y,x); // y=n+1-x
  mpz_set_ui(x,2); mpz_powm(x,x,y,n); // x=2^y mod n;

  cout << "trying phimat..." << x << endl;

  //mpz_add_ui(y,n,1); mpz_div_ui(y,y,2);
  unsigned long int z=0;
  while (z<1000000000)
    {
      int p=mpz_scan1(x,0); z+=p; 
      mpz_tdiv_q_2exp(x,x,p);
      //cout << z << ": " << x << endl;
      if (mpz_cmp_ui(x,1)==0) break;
      mpz_add(x,x,n);
    }

  cout << z << ":" << x << endl;
  
  if (mpz_cmp_ui(x,1)==0)
    {
      cout << "x=1 ... " << endl;
      mpz_mul_ui(x,n,4); mpz_sqrt(x,x); mpz_add_ui(x,x,z); // x=p+q
      mpz_mul(y,x,x);
      mpz_sub(y,y,n); mpz_sub(y,y,n); mpz_sub(y,y,n); mpz_sub(y,y,n);
      mpz_sqrt(y,y); // y=p+q
      mpz_add(y,x,y); mpz_div_ui(y,y,2); // y=p
      mpz_div(x,n,y); // x=q
      mpz_mul(x,x,y);
      cout << x << endl; 
      if (mpz_cmp(x,n)==0) // wirklich echte Teiler ermittelt?
        {
          mpz_div(x,n,y); // x=q
          if (mpz_probab_prime_p(x,probab_prime_checks))
            {
              cout << x << " is factor." << endl;
              Factorization_to_file << MAL(x) << flush;
            }
          else
            {
              cout << x << " is a composite factor." << endl;
              Factorization_to_file << MAL(x) << " [composite]" << flush;
            }
          if (mpz_probab_prime_p(y,probab_prime_checks))
            {
              cout << y << " is factor." << endl;
              Factorization_to_file << MAL(y) << flush;
            }
          else
            {
              cout << y << " is a composite factor." << endl;
              Factorization_to_file << MAL(y) << " [composite]" << flush;
            }      
          mpz_set_ui(n,1);
        }
    }
  
  mpz_clear(x); mpz_clear(y);
}

---