+
+/** Fibonacci number. The nth Fibonacci number F(n) is defined by the
+ * recurrence formula F(n)==F(n-1)+F(n-2) with F(0)==0 and F(1)==1.
+ *
+ * @param n an integer
+ * @return the nth Fibonacci number F(n) (an integer number)
+ * @exception range_error (argument must be an integer) */
+const numeric fibonacci(const numeric & n)
+{
+ if (!n.is_integer())
+ throw (std::range_error("numeric::fibonacci(): argument must be integer"));
+ // Method:
+ //
+ // This is based on an implementation that can be found in CLN's example
+ // directory. There, it is done recursively, which may be more elegant
+ // than our non-recursive implementation that has to resort to some bit-
+ // fiddling. This is, however, a matter of taste.
+ // The following addition formula holds:
+ //
+ // F(n+m) = F(m-1)*F(n) + F(m)*F(n+1) for m >= 1, n >= 0.
+ //
+ // (Proof: For fixed m, the LHS and the RHS satisfy the same recurrence
+ // w.r.t. n, and the initial values (n=0, n=1) agree. Hence all values
+ // agree.)
+ // Replace m by m+1:
+ // F(n+m+1) = F(m)*F(n) + F(m+1)*F(n+1) for m >= 0, n >= 0
+ // Now put in m = n, to get
+ // F(2n) = (F(n+1)-F(n))*F(n) + F(n)*F(n+1) = F(n)*(2*F(n+1) - F(n))
+ // F(2n+1) = F(n)^2 + F(n+1)^2
+ // hence
+ // F(2n+2) = F(n+1)*(2*F(n) + F(n+1))
+ if (n.is_zero())
+ return _num0();
+ if (n.is_negative())
+ if (n.is_even())
+ return -fibonacci(-n);
+ else
+ return fibonacci(-n);
+
+ ::cl_I u(0);
+ ::cl_I v(1);
+ ::cl_I m = The(::cl_I)(*n.value) >> 1L; // floor(n/2);
+ for (uintL bit=::integer_length(m); bit>0; --bit) {
+ // Since a squaring is cheaper than a multiplication, better use
+ // three squarings instead of one multiplication and two squarings.
+ ::cl_I u2 = ::square(u);
+ ::cl_I v2 = ::square(v);
+ if (::logbitp(bit-1, m)) {
+ v = ::square(u + v) - u2;
+ u = u2 + v2;
+ } else {
+ u = v2 - ::square(v - u);
+ v = u2 + v2;
+ }
+ }
+ if (n.is_even())
+ // Here we don't use the squaring formula because one multiplication
+ // is cheaper than two squarings.
+ return u * ((v << 1) - u);
+ else
+ return ::square(u) + ::square(v);
+}
+
+