* numeric::print(): increase precision for C source output.

[ginac.git] / ginac / numeric.cpp

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 55d5d556526e3038f47684e53a665527791cb9e8..0274d63af6fcee55c222c507253a7f124a53d8c5 100644 (file)
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -354,6 +354,11 @@ void numeric::print(const print_context & c, unsigned level) const
 
                std::ios::fmtflags oldflags = c.s.flags();
                c.s.setf(std::ios::scientific);
+               int oldprec = c.s.precision();
+               if (is_a<print_csrc_double>(c))
+                       c.s.precision(16);
+               else
+                       c.s.precision(7);
                if (this->is_rational() && !this->is_integer()) {
                        if (compare(_num0) > 0) {
                                c.s << "(";
@@ -381,6 +386,7 @@ void numeric::print(const print_context & c, unsigned level) const
                                c.s << to_double();
                }
                c.s.flags(oldflags);
+               c.s.precision(oldprec);
 
        } else {
                const std::string par_open  = is_a<print_latex>(c) ? "{(" : "(";
@@ -1546,14 +1552,13 @@ const numeric bernoulli(const numeric &nn)
 
        // 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];
+               return results[n/2-1];
 
-       results.reserve(n/2 + 1);
+       results.reserve(n/2);
        for (unsigned p=next_r; p<=n;  p+=2) {
                cln::cl_I  c = 1;  // seed for binonmial coefficients
                cln::cl_RA b = cln::cl_RA(1-p)/2;
@@ -1565,18 +1570,18 @@ const numeric bernoulli(const numeric &nn)
                if (p < (1UL<<cl_value_len/2)) {
                        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];
+                               b = b + c*results[k-1];
                        }
                } 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];
+                               b = b + c*results[k-1];
                        }
                }
                results.push_back(-b/(p+1));
        }
        next_r = n+2;
-       return results[n/2];
+       return results[n/2-1];
 }