-numeric bernoulli(numeric const & nn)
-{
- if (!nn.is_integer() || nn.is_negative())
- throw (std::range_error("numeric::bernoulli(): argument must be integer >= 0"));
- if (nn.is_zero())
- return numONE();
- if (!nn.compare(numONE()))
- return numeric(-1,2);
- if (nn.is_odd())
- return numZERO();
- // Until somebody has the Blues and comes up with a much better idea and
- // codes it (preferably in CLN) we make this a remembering function which
- // computes its results using the formula
- // B(nn) == - 1/(nn+1) * sum_{k=0}^{nn-1}(binomial(nn+1,k)*B(k))
- // whith B(0) == 1.
- static vector<numeric> results;
- static int highest_result = -1;
- int n = nn.sub(numTWO()).div(numTWO()).to_int();
- if (n <= highest_result)
- return results[n];
- if (results.capacity() < (unsigned)(n+1))
- results.reserve(n+1);
-
- numeric tmp; // used to store the sum
- for (int i=highest_result+1; i<=n; ++i) {
- // the first two elements:
- tmp = numeric(-2*i-1,2);
- // accumulate the remaining elements:
- for (int j=0; j<i; ++j)
- tmp += binomial(numeric(2*i+3),numeric(j*2+2))*results[j];
- // divide by -(nn+1) and store result:
- results.push_back(-tmp/numeric(2*i+3));
- }
- highest_result=n;
- return results[n];
+const numeric bernoulli(const numeric &nn)
+{
+ 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 formula below accomplishes this. It is a modification of the
+ // defining formula above but the computation of the binomial coefficients
+ // is carried along in an inline fashion. It also honors the fact that
+ // B_n is zero when n is odd and greater than 1.
+ //
+ // (There is an interesting relation with the tangent polynomials described
+ // in `Concrete Mathematics', which leads to a program a little faster 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.)
+
+ const unsigned n = nn.to_int();
+
+ // the special cases not covered by the algorithm below
+ if (n & 1)
+ return (n==1) ? _num_1_2 : _num0;
+ if (!n)
+ return _num1;
+
+ // store nonvanishing Bernoulli numbers here
+ static std::vector< cln::cl_RA > results;
+ static unsigned next_r = 0;
+
+ // algorithm not applicable to B(2), so just store it
+ if (!next_r) {
+ results.push_back(cln::recip(cln::cl_RA(6)));
+ next_r = 4;
+ }
+ if (n<next_r)
+ return results[n/2-1];
+
+ results.reserve(n/2);
+ for (unsigned p=next_r; p<=n; p+=2) {
+ cln::cl_I c = 1; // seed for binonmial coefficients
+ cln::cl_RA b = cln::cl_RA(1-p)/2;
+ const unsigned p3 = p+3;
+ const unsigned pm = p-2;
+ unsigned i, k, p_2;
+ // test if intermediate unsigned int can be represented by immediate
+ // objects by CLN (i.e. < 2^29 for 32 Bit machines, see <cln/object.h>)
+ if (p < (1UL<<cl_value_len/2)) {
+ for (i=2, k=1, p_2=p/2; i<=pm; i+=2, ++k, --p_2) {
+ c = cln::exquo(c * ((p3-i) * p_2), (i-1)*k);
+ b = b + c*results[k-1];
+ }
+ } else {
+ for (i=2, k=1, p_2=p/2; i<=pm; i+=2, ++k, --p_2) {
+ c = cln::exquo((c * (p3-i)) * p_2, cln::cl_I(i-1)*k);
+ b = b + c*results[k-1];
+ }
+ }
+ results.push_back(-b/(p+1));
+ }
+ next_r = n+2;
+ return results[n/2-1];
+}
+
+
+/** 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:
+ //
+ // 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);
+
+ cln::cl_I u(0);
+ cln::cl_I v(1);
+ cln::cl_I m = cln::the<cln::cl_I>(n.to_cl_N()) >> 1L; // floor(n/2);
+ for (uintL bit=cln::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.
+ cln::cl_I u2 = cln::square(u);
+ cln::cl_I v2 = cln::square(v);
+ if (cln::logbitp(bit-1, m)) {
+ v = cln::square(u + v) - u2;
+ u = u2 + v2;
+ } else {
+ u = v2 - cln::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 cln::square(u) + cln::square(v);