- constructor from strings once again accepts Lisp-style numbers like

[ginac.git] / ginac / numeric.cpp

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index c3ba81d6519719002be038f6d23441b593462321..f4dd738087de9fe0cc71b9f8a436694120b519e3 100644 (file)
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -231,7 +231,7 @@ numeric::numeric(const char *s) : basic(TINFO_numeric)
     // ss should represent a simple sum like 2+5*I
     std::string ss(s);
     // make it safe by adding explicit sign
-    if (ss.at(0) != '+' && ss.at(0) != '-')
+    if (ss.at(0) != '+' && ss.at(0) != '-' && ss.at(0) != '#')
         ss = '+' + ss;
     std::string::size_type delim;
     do {
@@ -1468,7 +1468,7 @@ const numeric binomial(const numeric & n, const numeric & k)
         }
     }
     
-    // should really be gamma(n+1)/(gamma(r+1)/gamma(n-r+1) or a suitable limit
+    // should really be gamma(n+1)/gamma(r+1)/gamma(n-r+1) or a suitable limit
     throw (std::range_error("numeric::binomial(): don´t know how to evaluate that."));
 }
 
@@ -1482,39 +1482,70 @@ const numeric bernoulli(const numeric & nn)
 {
     if (!nn.is_integer() || nn.is_negative())
         throw (std::range_error("numeric::bernoulli(): argument must be integer >= 0"));
-    if (nn.is_zero())
-        return _num1();
+    
+    // Method:
+    //
+    // The Bernoulli numbers are rational numbers that may be computed using
+    // the relation
+    //
+    //     B_n = - 1/(n+1) * sum_{k=0}^{n-1}(binomial(n+1,k)*B_k)
+    //
+    // with B(0) = 1.  Since the n'th Bernoulli number depends on all the
+    // previous ones, the computation is necessarily very expensive.  There are
+    // several other ways of computing them, a particularly good one being
+    // cl_I s = 1;
+    // cl_I c = n+1;
+    // cl_RA Bern = 0;
+    // for (unsigned i=0; i<n; i++) {
+    //     c = exquo(c*(i-n),(i+2));
+    //     Bern = Bern + c*s/(i+2);
+    //     s = s + expt_pos(cl_I(i+2),n);
+    // }
+    // return Bern;
+    // 
+    // 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().
+    // 
+    // (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
+    // addition to the remember table.  This doubles the memory footprint so
+    // we don't use it.)
+    
+    // the special cases not covered by the algorithm below
     if (!nn.compare(_num1()))
         return numeric(-1,2);
     if (nn.is_odd())
         return _num0();
-    // Until somebody has the blues and comes up with a much better idea and
-    // codes it (preferably in CLN) we make this a remembering function which
-    // computes its results using the defining formula
-    // B(nn) == - 1/(nn+1) * sum_{k=0}^{nn-1}(binomial(nn+1,k)*B(k))
-    // whith B(0) == 1.
-    // Be warned, though: the Bernoulli numbers are computationally very
-    // expensive anyhow and you shouldn't expect miracles to happen.
-    static vector<numeric> results;
-    static int highest_result = -1;
-    int n = nn.sub(_num2()).div(_num2()).to_int();
-    if (n <= highest_result)
-        return results[n];
-    if (results.capacity() < (unsigned)(n+1))
-        results.reserve(n+1);
     
-    numeric tmp;  // used to store the sum
-    for (int i=highest_result+1; i<=n; ++i) {
-        // the first two elements:
-        tmp = numeric(-2*i-1,2);
-        // accumulate the remaining elements:
-        for (int j=0; j<i; ++j)
-            tmp += binomial(numeric(2*i+3),numeric(j*2+2))*results[j];
-        // divide by -(nn+1) and store result:
-        results.push_back(-tmp/numeric(2*i+3));
+    // store nonvanishing Bernoulli numbers here
+    static vector< ::cl_RA > results;
+    static int highest_result = 0;
+    // algorithm not applicable to B(0), so just store it
+    if (results.size()==0)
+        results.push_back(::cl_RA(1));
+    
+    int n = nn.to_long();
+    for (int i=highest_result; i<n/2; ++i) {
+        ::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 = cl_I(n*m) * (B+results[j]) / (d1*d2);
+            n += 4;
+            m += 2;
+            d1 -= 1;
+            d2 -= 2;
+        }
+        B = (1 - ((B+1)/(2*i+3))) / (cl_I(1)<<(2*i+2));
+        results.push_back(B);
+        ++highest_result;
     }
-    highest_result=n;
-    return results[n];
+    return results[n/2];
 }
 
 
@@ -1528,8 +1559,16 @@ const numeric fibonacci(const numeric & n)
 {
     if (!n.is_integer())
         throw (std::range_error("numeric::fibonacci(): argument must be integer"));
+    // Method:
+    //
+    // This is based on an implementation that can be found in CLN's example
+    // directory.  There, it is done recursively, which may be more elegant
+    // than our non-recursive implementation that has to resort to some bit-
+    // fiddling.  This is, however, a matter of taste.
     // The following addition formula holds:
+    //
     //      F(n+m)   = F(m-1)*F(n) + F(m)*F(n+1)  for m >= 1, n >= 0.
+    //
     // (Proof: For fixed m, the LHS and the RHS satisfy the same recurrence
     // w.r.t. n, and the initial values (n=0, n=1) agree. Hence all values
     // agree.)