examples/fibonacci.cc

   1 // Compute and print the n-th Fibonacci number.
   2
   3 // We work with integers and real numbers.
   4 #include <cln/integer.h>
   5 #include <cln/real.h>
   6
   7 // We do I/O.
   8 #include <cln/io.h>
   9 #include <cln/integer_io.h>
  10
  11 // We use the timing functions.
  12 #include <cln/timing.h>
  13
  14 using namespace std;
  15 using namespace cln;
  16
  17 // F_n is defined through the recurrence relation
  18 //      F_0 = 0, F_1 = 1, F_(n+2) = F_(n+1) + F_n.
  19 // The following addition formula holds:
  20 //      F_(n+m)   = F_(m-1) * F_n + F_m * F_(n+1)  for m >= 1, n >= 0.
  21 // (Proof: For fixed m, the LHS and the RHS satisfy the same recurrence
  22 // w.r.t. n, and the initial values (n=0, n=1) agree. Hence all values agree.)
  23 // Replace m by m+1:
  24 //      F_(n+m+1) = F_m * F_n + F_(m+1) * F_(n+1)      for m >= 0, n >= 0
  25 // Now put in m = n, to get
  26 //      F_(2n) = (F_(n+1)-F_n) * F_n + F_n * F_(n+1) = F_n * (2*F_(n+1) - F_n)
  27 //      F_(2n+1) = F_n ^ 2 + F_(n+1) ^ 2
  28 // hence
  29 //      F_(2n+2) = F_(n+1) * (2*F_n + F_(n+1))
  30
  31 struct twofibs {
  32         cl_I u; // F_n
  33         cl_I v; // F_(n+1)
  34         // Constructor.
  35         twofibs (const cl_I& uu, const cl_I& vv) : u (uu), v (vv) {}
  36 };
  37
  38 // Returns F_n and F_(n+1). Assume n>=0.
  39 static const twofibs fibonacci2 (int n)
  40 {
  41         if (n==0)
  42                 return twofibs(0,1);
  43         int m = n/2; // floor(n/2)
  44         twofibs Fm = fibonacci2(m);
  45         // Since a squaring is cheaper than a multiplication, better use
  46         // three squarings instead of one multiplication and two squarings.
  47         cl_I u2 = square(Fm.u);
  48         cl_I v2 = square(Fm.v);
  49         if (n==2*m) {
  50                 // n = 2*m
  51                 cl_I uv2 = square(Fm.v - Fm.u);
  52                 return twofibs(v2 - uv2, u2 + v2);
  53         } else {
  54                 // n = 2*m+1
  55                 cl_I uv2 = square(Fm.u + Fm.v);
  56                 return twofibs(u2 + v2, uv2 - u2);
  57         }
  58 }
  59
  60 // Returns just F_n. Assume n>=0.
  61 const cl_I fibonacci (int n)
  62 {
  63         if (n==0)
  64                 return 0;
  65         int m = n/2; // floor(n/2)
  66         twofibs Fm = fibonacci2(m);
  67         if (n==2*m) {
  68                 // n = 2*m
  69                 // Here we don't use the squaring formula because
  70                 // one multiplication is cheaper than two squarings.
  71                 cl_I& u = Fm.u;
  72                 cl_I& v = Fm.v;
  73                 return u * ((v << 1) - u);
  74         } else {
  75                 // n = 2*m+1
  76                 cl_I u2 = square(Fm.u);
  77                 cl_I v2 = square(Fm.v);
  78                 return u2 + v2;
  79         }
  80 }
  81
  82 // The next routine is a variation of the above.  It is mathematically
  83 // equivalent but implemented in a non-recursive way.
  84 const cl_I fibonacci_compact (int n)
  85 {
  86         if (n==0)
  87                 return 0;
  88         cl_I u = 0;
  89         cl_I v = 1;
  90         cl_I m = n/2; // floor(n/2)
  91         for (uintL bit=integer_length(m); bit>0; --bit) {
  92                 // Since a squaring is cheaper than a multiplication, better use
  93                 // three squarings instead of one multiplication and two squarings.
  94                 cl_I u2 = square(u);
  95                 cl_I v2 = square(v);
  96                 if (logbitp(bit-1, m)) {
  97                         v = square(u + v) - u2;
  98                         u = u2 + v2;
  99                 } else {
 100                         u = v2 - square(v - u);
 101                         v = u2 + v2;
 102                 }
 103         }
 104         if (n==2*m)
 105                 // Here we don't use the squaring formula because
 106                 // one multiplication is cheaper than two squarings.
 107                 return u * ((v << 1) - u);
 108         else
 109                 return square(u) + square(v);
 110 }
 111
 112 // Returns just F_n, computed as the nearest integer to
 113 // ((1+sqrt(5))/2)^n/sqrt(5). Assume n>=0.
 114 const cl_I fibonacci_slow (int n)
 115 {
 116         // Need a precision of ((1+sqrt(5))/2)^-n.
 117         float_format_t prec = float_format((int)(0.208987641*n+5));
 118         cl_R sqrt5 = sqrt(cl_float(5,prec));
 119         cl_R phi = (1+sqrt5)/2;
 120         return round1( expt(phi,n)/sqrt5 );
 121 }
 122
 123 #ifndef TIMING
 124
 125 int main (int argc, char* argv[])
 126 {
 127         if (argc != 2) {
 128                 cerr << "Usage: fibonacci n" << endl;
 129                 return(1);
 130         }
 131         int n = atoi(argv[1]);
 132         cout << "fib(" << n << ") = " << fibonacci(n) << endl;
 133         return(0);
 134 }
 135
 136 #else // TIMING
 137
 138 int main (int argc, char* argv[])
 139 {
 140         int repetitions = 100;
 141         if ((argc >= 3) && !strcmp(argv[1],"-r")) {
 142                 repetitions = atoi(argv[2]);
 143                 argc -= 2; argv += 2;
 144         }
 145         if (argc != 2) {
 146                 cerr << "Usage: fibonacci n" << endl;
 147                 return(1);
 148         }
 149         int n = atoi(argv[1]);
 150         { CL_TIMING;
 151                 cout << "fib(" << n << ") = ";
 152                 for (int rep = repetitions-1; rep > 0; rep--)
 153                         fibonacci(n);
 154                 cout << fibonacci(n) << endl;
 155         }
 156         { CL_TIMING;
 157                 cout << "fib(" << n << ") = ";
 158                 for (int rep = repetitions-1; rep > 0; rep--)
 159                         fibonacci_compact(n);
 160                 cout << fibonacci_compact(n) << endl;
 161         }
 162         { CL_TIMING;
 163                 cout << "fib(" << n << ") = ";
 164                 for (int rep = repetitions-1; rep > 0; rep--)
 165                         fibonacci_slow(n);
 166                 cout << fibonacci_slow(n) << endl;
 167         }
 168         return(0);
 169 }
 170
 171 #endif