Make-safe patch to follow ... :)
Here it is, relative to the first patch and maybe a bit ugly. But using "(p < 32768)" does gain a speedup of about 5% for arguments smaller than approximately 1000. Cancelling 2 from both arguments of exquo changed the strange limiting value of the previous patch. Storing the binomial coefficients would certainly accelerate the calculation, again at the expense of something like doubled memory usage... Regards, Markus --- GiNaC-1.0.3_new_bernoulli/ginac/numeric.cpp Fri Jan 18 15:48:29 2002 +++ GiNaC-1.0.3_new_bernoulli_safe/ginac/numeric.cpp Mon Jan 21 18:02:03 2002 @@ -1534,24 +1534,35 @@ // algorithm not applicable to B(2), so just store it if (!next_r) { + results.push_back(); // results[0] is not used results.push_back(cln::recip(cln::cl_RA(6))); next_r = 4; } + if (n < next_r) + return results[n/2]; + + results.reserve(n/2 + 1); for (unsigned p = next_r; p <= n; p += 2) { - cln::cl_I c = 1; + cln::cl_I c = 1; // binonmial coefficients 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]; - } + unsigned i, k, p_2; + // test if intermediate unsigned int results < 2^29 + if (p < 32768) + 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]; + } + 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]; + } results.push_back(- b / (p + 1)); - next_r += 2; - } - return results[n/2 - 1]; + } + next_r = n + 2; + return results[n/2]; }