-numeric doublefactorial(numeric const & nn)
-{
- // We store the results separately for even and odd arguments. This has
- // the advantage that we don't have to compute any even result at all if
- // the function is always called with odd arguments and vice versa. There
- // is no tradeoff involved in this, it is guaranteed to save time as well
- // as memory. (If this is not enough justification consider the Gamma
- // function of half integer arguments: it only needs odd doublefactorials.)
- static vector<numeric> evenresults;
- static int highest_evenresult = -1;
- static vector<numeric> oddresults;
- static int highest_oddresult = -1;
-
- if ( nn == numeric(-1) ) {
- return numONE();
- }
- if ( !nn.is_nonneg_integer() ) {
- throw (std::range_error("numeric::doublefactorial(): argument must be integer >= -1"));
- }
- if ( nn.is_even() ) {
- int n = nn.div(numTWO()).to_int();
- if ( n <= highest_evenresult ) {
- return evenresults[n];
- }
- if ( evenresults.capacity() < (unsigned)(n+1) ) {
- evenresults.reserve(n+1);
- }
- if ( highest_evenresult < 0 ) {
- evenresults.push_back(numONE());
- highest_evenresult=0;
- }
- for (int i=highest_evenresult+1; i<=n; i++) {
- evenresults.push_back(numeric(evenresults[i-1].mul(numeric(i*2))));
- }
- highest_evenresult=n;
- return evenresults[n];
- } else {
- int n = nn.sub(numONE()).div(numTWO()).to_int();
- if ( n <= highest_oddresult ) {
- return oddresults[n];
- }
- if ( oddresults.capacity() < (unsigned)n ) {
- oddresults.reserve(n+1);
- }
- if ( highest_oddresult < 0 ) {
- oddresults.push_back(numONE());
- highest_oddresult=0;
- }
- for (int i=highest_oddresult+1; i<=n; i++) {
- oddresults.push_back(numeric(oddresults[i-1].mul(numeric(i*2+1))));
- }
- highest_oddresult=n;
- return oddresults[n];
- }
-}
-
-/** The Binomial function. It computes the binomial coefficients. If the
- * arguments are both nonnegative integers and 0 <= k <= n, then
- * binomial(n, k) = n!/k!/(n-k)! which is the number of ways of choosing k
- * objects from n distinct objects. If k > n, then binomial(n,k) returns 0. */
-numeric binomial(numeric const & n, numeric const & k)
-{
- if (n.is_nonneg_integer() && k.is_nonneg_integer()) {
- return numeric(binomial(n.to_int(),k.to_int())); // -> CLN
- } else {
- // should really be gamma(n+1)/(gamma(r+1)/gamma(n-r+1)
- return numeric(0);
- }
- // return factorial(n).div(factorial(k).mul(factorial(n.sub(k))));
+const numeric doublefactorial(const numeric &n)
+{
+ if (n.is_equal(_num_1))
+ return _num1;
+
+ if (!n.is_nonneg_integer())
+ throw std::range_error("numeric::doublefactorial(): argument must be integer >= -1");
+
+ return numeric(cln::doublefactorial(n.to_int()));
+}
+
+
+/** The Binomial coefficients. It computes the binomial coefficients. For
+ * integer n and k and positive n this is the number of ways of choosing k
+ * objects from n distinct objects. If n is negative, the formula
+ * binomial(n,k) == (-1)^k*binomial(k-n-1,k) is used to compute the result. */
+const numeric binomial(const numeric &n, const numeric &k)
+{
+ if (n.is_integer() && k.is_integer()) {
+ if (n.is_nonneg_integer()) {
+ if (k.compare(n)!=1 && k.compare(_num0)!=-1)
+ return numeric(cln::binomial(n.to_int(),k.to_int()));
+ else
+ return _num0;
+ } else {
+ return _num_1.power(k)*binomial(k-n-_num1,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) */
+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);