Alternative bernoulli(const numeric &)

[Markus Nullmeier]

> > 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];
 }