+ if (n.is_integer() && k.is_integer()) {
+ if (n.is_nonneg_integer()) {
+ if (k.compare(n)!=1 && k.compare(numZERO())!=-1)
+ return numeric(::binomial(n.to_int(),k.to_int())); // -> CLN
+ else
+ return numZERO();
+ } else {
+ return numMINUSONE().power(k)*binomial(k-n-numONE(),k);
+ }
+ }
+
+ // should really be gamma(n+1)/(gamma(r+1)/gamma(n-r+1) or a suitable limit
+ throw (std::range_error("numeric::binomial(): donĀ“t know how to evaluate that."));
+}
+
+/** Bernoulli number. The nth Bernoulli number is the coefficient of x^n/n!
+ * in the expansion of the function x/(e^x-1).
+ *
+ * @return the nth Bernoulli number (a rational number).
+ * @exception range_error (argument must be integer >= 0) */
+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));