+ if (!nn.is_integer() || nn.is_negative())
+ throw (std::range_error("numeric::bernoulli(): argument must be integer >= 0"));
+
+ // Method:
+ //
+ // The Bernoulli numbers are rational numbers that may be computed using
+ // the relation
+ //
+ // B_n = - 1/(n+1) * sum_{k=0}^{n-1}(binomial(n+1,k)*B_k)
+ //
+ // with B(0) = 1. Since the n'th Bernoulli number depends on all the
+ // previous ones, the computation is necessarily very expensive. There are
+ // several other ways of computing them, a particularly good one being
+ // cl_I s = 1;
+ // cl_I c = n+1;
+ // cl_RA Bern = 0;
+ // for (unsigned i=0; i<n; i++) {
+ // c = exquo(c*(i-n),(i+2));
+ // Bern = Bern + c*s/(i+2);
+ // s = s + expt_pos(cl_I(i+2),n);
+ // }
+ // return Bern;
+ //
+ // But if somebody works with the n'th Bernoulli number she is likely to
+ // also need all previous Bernoulli numbers. So we need a complete remember
+ // table and above divide and conquer algorithm is not suited to build one
+ // up. The code below is adapted from Pari's function bernvec().
+ //
+ // (There is an interesting relation with the tangent polynomials described
+ // in `Concrete Mathematics', which leads to a program twice as fast as our
+ // implementation below, but it requires storing one such polynomial in
+ // addition to the remember table. This doubles the memory footprint so
+ // we don't use it.)
+
+ // the special cases not covered by the algorithm below
+ if (!nn.compare(_num1()))
+ return numeric(-1,2);
+ if (nn.is_odd())
+ return _num0();
+
+ // store nonvanishing Bernoulli numbers here
+ static vector< ::cl_RA > results;
+ static int highest_result = 0;
+ // algorithm not applicable to B(0), so just store it
+ if (results.size()==0)
+ results.push_back(::cl_RA(1));
+
+ int n = nn.to_long();
+ for (int i=highest_result; i<n/2; ++i) {
+ ::cl_RA B = 0;
+ long n = 8;
+ long m = 5;
+ long d1 = i;
+ long d2 = 2*i-1;
+ for (int j=i; j>0; --j) {
+ B = cl_I(n*m) * (B+results[j]) / (d1*d2);
+ n += 4;
+ m += 2;
+ d1 -= 1;
+ d2 -= 2;
+ }
+ B = (1 - ((B+1)/(2*i+3))) / (cl_I(1)<<(2*i+2));
+ results.push_back(B);
+ ++highest_result;
+ }
+ return results[n/2];
+}
+
+
+/** 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;
+ }