Hello, being unware of GiNaC's implemetation, I wrote my own Bernoulli number function after reading Peter Stone's straightforward description in his Maple Worksheet http://minyos.its.rmit.edu.au/~e05254/MAPLE/ENGMATH/TAYLOR/BERNUM.MWS Later I adapted my function to bernoulli() in numeric.cpp, as it turned out to be about twice as fast as GiNaC 1.0.3 (gcc 2.95, -O2). Thus maybe you'll find the enclosed patch useful in some way or the other. Hopefully no hidden bugs remain. Markus --- GiNaC-1.0.3/ginac/numeric.cpp Wed Dec 19 12:11:46 2001 +++ GiNaC-1.0.3_new_bernoulli/ginac/numeric.cpp Fri Jan 18 14:44:16 2002 @@ -1512,7 +1512,7 @@ // 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. // // (There is an interesting relation with the tangent polynomials described // in `Concrete Mathematics', which leads to a program twice as fast as our @@ -1520,38 +1520,38 @@ // addition to the remember table. This doubles the memory footprint so // we don't use it.) + 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; i<n/2; ++i) { - cln::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 = cln::cl_I(n*m) * (B+results[j]) / (d1*d2); - n += 4; - m += 2; - d1 -= 1; - d2 -= 2; - } - B = (1 - ((B+1)/(2*i+3))) / (cln::cl_I(1)<<(2*i+2)); - results.push_back(B); - ++highest_result; - } - return results[n/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; + unsigned p3 = p+3; + unsigned p2 = p+2; + 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]; + } + results.push_back(- b / (p + 1)); + next_r += 2; + } + return results[n/2 - 1]; }