- status: open --> closed
The zeta function provides incorrect numerical output for numbers that aren't even integers greater than some number between 33.44 and 33.45. The output for 33.44 is 1.000000000128721, and for 33.45 is 0.0. It should be clear from the definition ( http://dlmf.nist.gov/25.2 ) that for finite positive real values zeta should be > 1 (>= 1.0 given the accuracy limitations of floating point). For values this large, direct calculation using the definition should suffice.
Another potentially useful thing to do would be to provide a function zetam1(s) that is formally equal to zeta(s) - 1 that can work for large values of s where the difference between zeta(s) and 1 would be lost in floating point cancellation.
I've just fixed this bug and pushed the fix to master. The mistake was that we were incorrectly calculating the number of terms to use in our series approximation.
Re the zetam1 proposal, this should probably be discussed on the mailing list. Note that, using bigfloats, we can represent floating point numbers arbitrarily close to 1. As such, I think a "zetam1" function would only be useful if it's much cheaper to compute zeta(s)-1 than to compute zeta(s). And of course that's not going to be the case. Other people's opions might differ... :-)
I think zetam1 is useful, but that's at least a feature request, not a bug.
If you're computing zeta - 1, a zetam1 function could allow you do the work with doubles instead of bigfloats, or less precise bigfloats than would otherwise would be needed to compute zeta - 1 to decent accuracy.
Maybe I've misunderstood, but doesn't that only help if you've got a way to calculate zeta(s)-1? I mean, yes, you can calculate zeta(s) with higher floating precision then subtract 1 and return the result in just a double, but that seems pretty niche...
I haven't seen any papers that say "Ahah, near the real line, you can compute zeta(x)-1 with the following series" (and then don't do something silly like adding and subtracting terms of order one), but maybe I've not been looking in the right place?
I'm sure there are better ways, but you can do sum(1/n^s, n, 1, inf) via a Euler-Maclaurin series to get zetam1(s). The Euler-Maclaurin series converges quite a bit faster than the original series.
Great! So on what region should we be using that method? Using it, you can calculate zeta by just adding one at the end. (which is why I was/am slightly sceptical)
I note that the paper that describes the series we're using starts by talking about Euler-Maclaurin summation of the standard series and explains that their method is faster. But maybe this is something about real vs. complex arguments?
I believe Euler-Maclaurin series is valid over the entire plane, but I'd have to go back read over everything again.
But as I mentioned zetam1 is a feature request.