From: Richard Kreckel Date: Fri, 18 Jan 2002 19:59:09 +0000 (+0000) Subject: * Faster Bernoulli numbers (Markus Nullmeier). X-Git-Tag: release_1-0-4~7 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=c86cff51ac5a42f86387ef7bc767f1274137350b * Faster Bernoulli numbers (Markus Nullmeier). --- diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 01c58453..e621607a 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -1488,7 +1488,7 @@ 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 @@ -1512,46 +1512,49 @@ const numeric bernoulli(const numeric &nn) // 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(). + // 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 twice as fast as our - // implementation below, but it requires storing one such polynomial in + // 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 (nn.is_equal(_num1)) - return _num_1_2; - if (nn.is_odd()) - return _num0; - + 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 int highest_result = 0; - // algorithm not applicable to B(0), so just store it - if (results.empty()) - results.push_back(cln::cl_RA(1)); - - int n = nn.to_long(); - for (int i=highest_result; i0; --j) { - B = cln::cl_I(n*m) * (B+results[j]) / (d1*d2); - n += 4; - m += 2; - d1 -= 1; - d2 -= 2; + 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; + } + for (unsigned p=next_r; p<=n; p+=2) { + cln::cl_I c = 1; + cln::cl_RA b = cln::cl_RA(1-p)/2; + const unsigned p3 = p+3; + const unsigned p2 = p+2; + const unsigned pm = p-2; + unsigned i, k; + for (i=2, k=0; i <= pm; i += 2, k++) { + c = cln::exquo(c * ((p3 - i)*(p2 - i)), (i - 1)*i); + b = b + c * results[k]; } - B = (1 - ((B+1)/(2*i+3))) / (cln::cl_I(1)<<(2*i+2)); - results.push_back(B); - ++highest_result; + results.push_back(-b / (p+1)); + next_r += 2; } - return results[n/2]; + return results[n/2 - 1]; }