Alternative bernoulli(const numeric &)

[Markus Nullmeier]

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