## #2682 Function zeta fails numerically for large numbers that aren't even integers

None
closed
nobody
None
5
2014-01-28
2014-01-23
Sean Lake
No

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.

## Discussion

• Rupert Swarbrick - 2014-01-26
• status: open --> closed

• Rupert Swarbrick - 2014-01-26

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... :-)

• Raymond Toy - 2014-01-26

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.

• Rupert Swarbrick - 2014-01-26

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?

• Raymond Toy - 2014-01-28

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.

• Rupert Swarbrick - 2014-01-28

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?

• Raymond Toy - 2014-01-28

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.