1 // Compute and print the n-th Fibonacci number.
3 // We work with integers and real numbers.
4 #include <cln/integer.h>
9 #include <cln/integer_io.h>
11 // We use the timing functions.
12 #include <cln/timing.h>
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.)
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
29 // F_(2n+2) = F_(n+1) * (2*F_n + F_(n+1))
35 twofibs (const cl_I& uu, const cl_I& vv) : u (uu), v (vv) {}
38 // Returns F_n and F_(n+1). Assume n>=0.
39 static const twofibs fibonacci2 (int n)
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);
51 cl_I uv2 = square(Fm.v - Fm.u);
52 return twofibs(v2 - uv2, u2 + v2);
55 cl_I uv2 = square(Fm.u + Fm.v);
56 return twofibs(u2 + v2, uv2 - u2);
60 // Returns just F_n. Assume n>=0.
61 const cl_I fibonacci (int n)
65 int m = n/2; // floor(n/2)
66 twofibs Fm = fibonacci2(m);
69 // Here we don't use the squaring formula because
70 // one multiplication is cheaper than two squarings.
73 return u * ((v << 1) - u);
76 cl_I u2 = square(Fm.u);
77 cl_I v2 = square(Fm.v);
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)
90 cl_I m = n/2; // floor(n/2)
91 for (uintC 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.
96 if (logbitp(bit-1, m)) {
97 v = square(u + v) - u2;
100 u = v2 - square(v - u);
105 // Here we don't use the squaring formula because
106 // one multiplication is cheaper than two squarings.
107 return u * ((v << 1) - u);
109 return square(u) + square(v);
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)
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 );
125 int main (int argc, char* argv[])
128 cerr << "Usage: fibonacci n" << endl;
131 int n = atoi(argv[1]);
132 cout << "fib(" << n << ") = " << fibonacci(n) << endl;
138 int main (int argc, char* argv[])
140 int repetitions = 100;
141 if ((argc >= 3) && !strcmp(argv[1],"-r")) {
142 repetitions = atoi(argv[2]);
143 argc -= 2; argv += 2;
146 cerr << "Usage: fibonacci n" << endl;
149 int n = atoi(argv[1]);
151 cout << "fib(" << n << ") = ";
152 for (int rep = repetitions-1; rep > 0; rep--)
154 cout << fibonacci(n) << endl;
157 cout << "fib(" << n << ") = ";
158 for (int rep = repetitions-1; rep > 0; rep--)
159 fibonacci_compact(n);
160 cout << fibonacci_compact(n) << endl;
163 cout << "fib(" << n << ") = ";
164 for (int rep = repetitions-1; rep > 0; rep--)
166 cout << fibonacci_slow(n) << endl;