fibonacci_ppm1.cc

Go to the documentation of this file.

#include "polphi_template.H"

class Tfibpair
{
 private:
  mpz_t x1,x2; // (x1,x2) stands for (x1*lambda + x2), which complies with (f[n],f[n-1])
  mpz_t h1,h2,h3;

  inline const Tfibpair& operator= (const Tfibpair&) const { MARK; exit(9); return *this; }
    // dummy to forbid any implicit inherited use of the assignment operator

 public:
  Tfibpair()
   {
     // Default-Startwert ist lambda^0 = 1 = 0*lambda+1 = (0,1), also f1=F[0]=0, f0=F[-1]=1
     mpz_init_set_ui(x1,0);
     mpz_init_set_ui(x2,1);
     mpz_init(h1); mpz_init(h2); mpz_init(h3);
   }
  Tfibpair(const unsigned int f1, const unsigned int f0)
   {
     // Default-Startwert ist lambda^0 = 1 = 0*lambda+1 = (0,1), also f1=F[0]=0, f0=F[-1]=1
     // Für den Startwert der Fibonaccifolge wähle man f1=F[1]=1, f0=F[0]=0
     mpz_init_set_ui(x1,f1);
     mpz_init_set_ui(x2,f0);
     mpz_init(h1); mpz_init(h2); mpz_init(h3);
   }
  Tfibpair(const Tfibpair &X)
   {
     mpz_init_set(x1,X.x1);
     mpz_init_set(x2,X.x2);
     mpz_init(h1); mpz_init(h2); mpz_init(h3);
   }
  ~Tfibpair()
   {
     mpz_clear(h3); mpz_clear(h2); mpz_clear(h1);
     mpz_clear(x2); mpz_clear(x1);
   }

  inline void set(const unsigned int f1, const unsigned int f0)
   {
     mpz_set_ui(x1,f1);
     mpz_set_ui(x2,f0);
   }
  inline void set(const signed int f1, const signed int f0)
   {
     mpz_set_si(x1,f1);
     mpz_set_si(x2,f0);
   }
  inline void set(const Tfibpair &X)
   {
     mpz_set(x1,X.x1);
     mpz_set(x2,X.x2);
     // Hilfsvariablen bleiben unberücksichtigt!
   }


  inline const mpz_t& fn(void) const { return x1; }
  inline const mpz_t& fnm1(void) const { return x2; }
  inline const mpz_t& Ln(void) { mpz_add(h1,x2,x2); mpz_add(h1,h1,x1); return h1; }

  inline void mod(const mpz_t m)
   {
     mpz_mod(x1,x1,m);
     mpz_mod(x2,x2,m);
   }
  inline void mul(const Tfibpair &Q)
   {
     mpz_add(h1,x1,x2); mpz_add(h2,Q.x1,Q.x2); mpz_mul(h1,h1,h2);
     mpz_mul(h2,x1,Q.x1); 
     mpz_mul(h3,x2,Q.x2); 
     mpz_sub(x1,h1,h3); mpz_add(x2,h2,h3);
   }
  inline void square() 
   {
     mpz_add(h1,x1,x2); mpz_mul(h1,h1,h1);
     mpz_mul(h2,x1,x1); 
     mpz_mul(h3,x2,x2); 
     mpz_sub(x1,h1,h3); mpz_add(x2,h2,h3);
   }
  inline void fastsquare()
   {
     // assumes that first index is already even
#ifdef DEBUG
     const bool flag = (mpz_cmp_ui(x1,1)==0 && mpz_cmp_ui(x2,0)==0);
     if (flag) cout << "error, error!" << endl;
#endif
     mpz_mul(h1,x1,x1); mpz_mul(h2,x2,x2);
     mpz_add(x2,h1,h2);
     mpz_mul_2exp(h1,h1,2); mpz_sub(h1,h1,h2); mpz_add_ui(h1,h1,2);
     mpz_sub(x1,h1,x2);
   }
  inline void pow(const unsigned int n)
   {
     // calculating (x1*lambda+x2)^n
     Tfibpair Q(*this);
     if ((n&1)==0) set(0,1);
     for (unsigned int i=2; i<=n; i<<=1)
      {
        Q.square();
        if (n&i) mul(Q);
      }
   }
  inline void powmod(const unsigned int n, const mpz_t m)
   {
     // standard binary powering algorithm
     // this one goes from the lowest bit upwards
     Tfibpair Q(*this);
     if ((n&1)==0) set(0,1);
     for (unsigned int i=2; i<=n; i<<=1)
      {
        Q.square(); Q.mod(m);
        if (n&i) { mul(Q); mod(m); }
      }
   }

#if 0
  inline void fastpowmod(const unsigned int n, const mpz_t m)
   {
     // binary powering algorithm, just the other way round:
     // this one goes from highest bit downwards.

     // assumes, that
     // (1) index of the existing fibonacci number is even, so that
     //     assumptions on Lucasvalues are correct
     // (2) n is a positive value, so that at least one bit is set

     const Tfibpair saved(*this);
     signed int b = 8*sizeof(b)-1; // highest possible bitposition
     while ( (n&(1<<b))==0 ) --b; // find highest set bit
     while (--b>=0)
      {
        fastsquare(); mod(m);
        if (n&(1<<b)) { mul(saved); mod(m); }
      }
   }
#else
  inline void fastpowmod(const unsigned int n, const mpz_t m)
   {
     // binary powering algorithm, just the other way round:
     // this one goes from highest bit downwards.

     // assumes, that
     // (1) index of the existing fibonacci number is even, so that
     //     assumptions on Lucasvalues are correct
     // (2) n is a positive value, so that at least one bit is set

     mpz_t dx1, dx2;
     mpz_init_set(dx1,x1); mpz_init_set(dx2,x2); mpz_add(h3,x1,x2);
     signed int b = 8*sizeof(b)-1; // highest possible bitposition
     while ( (n&(1<<b))==0 ) --b; // find highest set bit
     while (--b>=0)
      {
        //fastsquare();
        mpz_mul(h1,x1,x1); mpz_mul(h2,x2,x2);
        mpz_add(x2,h1,h2);
        mpz_mul_2exp(h1,h1,2); mpz_sub(h1,h1,h2); mpz_add_ui(h1,h1,2);
        mpz_sub(x1,h1,x2);
        mod(m);
        if (n&(1<<b))
         {
           //mul(saved);
           //mpz_add(h1,x1,x2); already done!
           //mpz_add(h2,dx1,dx2); // in h3!
           mpz_mul(h1,h1,h3);
           mpz_mul(h2,x2,dx2);
           mpz_mul(x2,x1,dx1);  
           mpz_sub(x1,h1,h2); mpz_add(x2,x2,h2);
           mod(m);
         }
      }
     mpz_clear(dx2); mpz_clear(dx1);    
   }
#endif

#if 1
  // the following class functions are untested, but they should work:

  inline void set_Fibonacci(const signed int n)
   {
     if (n==0) { set(0,1); return; }
     if (n<0) { set(1,-1); pow(-n); } else { set(1,0); pow(n); }
   }

  inline void step_forward(unsigned int d)
   {
     // (F[n],F[n-1]) -> (F[n+d],F[n+d-1])
     // one can also use (F[n+d],F[n+d-1])=(F[n],F[n-1])*(F[d],F[d-1])
     // (1) check threshold, whether powering will be worthwhile
     // (2) compute (F[d],F[d-1])=lambda^d and multiply...
     if (d<100)
      {
        while(d--) { mpz_swap(x1,x2); mpz_add(x1,x1,x2); }
      }
     else
      {
        Tfibpair temp(1,0);
        temp.pow(d);
        mul(temp); // "bigstep d" ;-)
      }
   }
  inline void step_backward(unsigned int d)
   {
     // (F[n],F[n-1]) -> (F[n-d],F[n-d-1])
     // one can also use (F[n+d],F[n+d-1])=(F[n],F[n-1])*(F[d],F[d-1])
     // (1) check threshold, whether powering will be worthwhile
     // (2) compute (F[-d],F[-d-1])=lambda^(-d) and multiply...
     if (d<100)
      {
        while(d--) { mpz_sub(x1,x1,x2); mpz_swap(x1,x2); }
      }
     else
      {
        Tfibpair temp;
        temp.set_Fibonacci(-1); 
        temp.pow(d);
        mul(temp); // "bigstep -d" ;-)
      }
   }
#endif

};


class CRingFib
{
private:
 Tfibpair a;

  inline const CRingFib& operator= (const CRingFib&) const { MARK; exit(9); return *this; }
    // dummy to forbid any implicit inherited use of the assignment operator

public:
 static const std::string name;
 static const int SieveIntervalSize=2*3*5*7*11*13*17;
 CRingFib() : a() { }
 CRingFib(const CRingFib &x) : a(x.a) { }
 ~CRingFib() { }
 bool set_startvalue(const unsigned int seed=1)
  {
    a.set(1,0);
    return seed==1; // reject seed!=1
  }
 void set(const CRingFib &x) { a.set(x.a); }
 void pow_mod(const unsigned int i, const mpz_t M) { a.powmod(i,M); }
 void fast_pow_mod(const unsigned int i, const mpz_t M) { a.fastpowmod(i,M); } // see above
 void test_gcd(mpz_t x, const mpz_t M) { mpz_gcd(x,a.fn(),M); }
 friend class CRingFibPhase2;
};
const std::string CRingFib::name = "fib";


class CRingFibPhase2 : private ForbidAssignment
{
private:
 Tfibpair F, stepF, mulF;
 mpz_t Fii, Lii, Liip, FD, LD, temp1, temp2;
public:
 explicit CRingFibPhase2(const CRingFib &_F) : F(_F.a), stepF(), mulF()
  {
    F.powmod(2,n); // um wirklich sicherzustellen, dass Index gerade ist...
    stepF.set(F), mulF.set(F);
    mpz_init(Fii); mpz_init(Lii); mpz_init(Liip);
    mpz_init(FD); mpz_init(LD); mpz_init(temp1); mpz_init(temp2);
  }
 ~CRingFibPhase2()
  {
    mpz_clear(Fii); mpz_clear(Lii); mpz_clear(Liip);
    mpz_clear(FD); mpz_clear(LD); mpz_clear(temp1); mpz_clear(temp2);
  }

 void get_polynomdef_point(mpz_t x)
  {
    mpz_set(x,stepF.fn()); mpz_mul(x,x,x); mpz_mod(x,x,n);
  }
 void calc_polynomdef_next_point()
  {
    stepF.mul(mulF); stepF.mod(n);
  }

 void calc_EvalStartingPoint(const int D, const double ii)
  {
    Tfibpair h(F); // only needed temporarily
    h.powmod(D,n);
    mpz_set(FD,h.fn()); // Fibonacciwert F[D]
    mpz_set(LD,h.Ln()); // Lucaswert L[D]

    Tfibpair hp(h);
    int pow = static_cast<int>(ii/D);
    h.powmod(pow,n);
    mpz_set(Fii,h.fn());
    mpz_set(Lii,h.Ln());
 
    if (--pow<0) pow=-pow;
    hp.powmod(pow,n);
    mpz_set(Liip,hp.Ln()); // Lucaswert L[ii-D]
  }

 void get_point_and_calc_next_point(mpz_t x)
  {
    mpz_mul(x,Fii,Fii); mpz_mod(x,x,n);
    // compute values for the next interval...
    //ii+=D; // neues Intervall (bis i haben wir nun abgearbeitet)

    // F[ii+D]=(F[ii]*L[D]+F[D]*L[ii])/2
    // L[ii+D]=L[ii]*L[D]-L[ii-D]

    mpz_mul(temp1,Fii,LD); mpz_mul(temp2,Lii,FD);
    mpz_add(Fii,temp1,temp2); mpz_mod(Fii,Fii,n);
    if (mpz_odd_p(Fii)) mpz_add(Fii,Fii,n);
    mpz_tdiv_q_2exp(Fii,Fii,1);

    mpz_mul(temp1,Lii,LD); mpz_sub(temp1,temp1,Liip); mpz_mod(Liip,temp1,n);
    mpz_swap(Liip,Lii);
  }
  
};



// explicit instantiation (not necessary)
template void polphi_template<CRingFib,CRingFibPhase2>(const int phase1, const double phase2);

#define fibonacci_ppm1(phase1,phase2) polphi_template<CRingFib,CRingFibPhase2>(phase1,phase2)

---