[#125] Implement Hurwitz zeta function

[rockbmb]

Closes #125. Still a WIP.

[rockbmb]

@Bodigrim is it worth it to implement the Hurwitz zeta function with both real arguments? It'll take some more work as it'll require the gamma function, which I don't think is implemented.

The Riemann zeta and Dirichlet beta will only need an integral first argument, so for a first step that will be sufficient.

[Bodigrim]

Gamma function is available in gamma package, but only for Double and Float. It would be better to have it for any Floating.

I think as a first stab we may keep the first argument integer.

[rockbmb]

but only for Double and Float.

Yikes, that won't do. I'll look into submitting a PR to change that, although I don't know if it'll happen. For now, as you agreed, first integer is Integral.

[rockbmb]

@Bodigrim I refactored the code to instead use a simple implementation of the ideas mentioned in the thesis (formula 4.8.9 of linked page) I mentioned while working on #120.

The Hurwitz zeta function is now stabler, but there are some issues with the precision tests for the Riemann zeta and Dirichlet beta tests. Will delete them if they no longer make sense for Hurwitz zeta.

[Bodigrim]

Let us name it zetaHurwitz or something else, different from zeta. Just to avoid confusion with Riemann zeta.

[Bodigrim]

zetasProperty2 is tough to satisfy (I remember having extremely hard times with it), but it is the only way to check the precision, which zetas was asked for.

[rockbmb]

@Bodigrim I think the reason the tests fail is that with 4.8.5 is because in the "Evaluating the error bound" section, it is stated that

For a precision of P bits, we should choose N ∼ M ∼ P

However, the current method (disregard the comments, outdated. Will fix.) to set N and M is crude, just relying on the fact that P bits of precision will result in P * (log 2 / log 10) ~ P * 3.3 decimal digits of precision. I keep getting a lot of NaNs, which may be a sign some (truncated) sums are consuming more terms than necessary.

I think my next step will be to implement the error bounding described in Theorem 4.8.1, and from this bounding derive better N, M. Ultimately this method won't be perfect either because we don't have the "ball" arbitrary precision types that the Arb library does.

[Bodigrim]

It looks to me that tests fail for eps close or smaller then 2^-52. In such case it is just a numeric noise, coming from Double itself, which cannot be cured by computing more terms of a series or whatever. I think it makes sense to restrict eps >= 2^-50 in these tests.

[rockbmb]

@Bodigrim You're right. It's not worth asking for precision greater than 1e-52 on a Double, or 1e-50 for Fixed Prec50, for example.

However, I'll still look into this issue with error bounding, and also the way and order the sums are calculated in. The previous implementation also received such high precisions in this test (i.e. 2e-67) and passed without generating improper values, so hitting NaN so quickly is a red flag for this implementation.

[rockbmb]

@Bodigrim now that the Hurwitz zeta function is working for integer arguments, I'm going to try extending it to real arguments, first and then complex. This requires the gamma function, which has been implemented in the package we mentioned before. It seems abandoned, so first I'll just write it locally and then consider a PR.

[Bodigrim]

Good!

I looked at gamma package more closely. Good news: it still compiles. Bad news: not only gamma itself looks unsupported, but it also depends on two more presumably unsupported packages.

[rockbmb]

@Bodigrim there are other issues, such as the relative complexity of the implementation and it's lack of documentation. I'm sure there are a lot of very clever tricks used throughout, but without an explanation it's hard to maintain them.

[rockbmb]

@Bodigrim Apologies, I made an incorrect statement: it is not necessary to use the gamma function for an extended version of zetaHurwitz. The Pochhammer symbol (for the rising factorial) has the same definition for either integer or real values.

However, turning the integer argument into a floating argument reduces precision by default, simple tests like beta(3) and zeta(2) fail. Not sure why that is, could have to do with numerical stability of Prelude functions like (**) or maybe a previous assumption I made in the code that is no longer true... This PR is turning into quite the task.

There's one more matter I'd like to discuss. While the gamma package will not be useful to us this time, I'm certain it will be in the future because of its importance in analytical number theory. Do you find it worthwhile to bring the package up to speed, for future utilization? At least an initial PR, to ask the maintainer for his opinion.

[Bodigrim]

Yeah, dealing with numerical stability is tough.

---

With regards to gamma package. From my perspective it would be nice to migrate https://github.com/mokus0/gamma/blob/master/extras/LanczosConstants.hs to any existing matrix package (actually, matrix will do) and release it as a user-facing module of gamma.

As a zero step, one can add gamma and its dependencies to Stackage.

I typically write to maintainers privately, asking if they have any plans of future development, explaining my intention to contribute and to help with life cycling. People usually like when someone finds their abandoned code useful and is willing to continue development.

[rockbmb]

@Bodigrim

migrate https://github.com/mokus0/gamma/blob/master/extras/LanczosConstants.hs to any existing matrix package
and release it as a user-facing module of gamma

I'm not sure I understand the idea here. Can you explain what your plan for the LanczosConstants module would be?

[Bodigrim]

AFAIU from skimming through gamma, it uses Lancsoz approximation to compute values of the gamma function. Namely, there are two (one for Float and one for Double) precomputed lists of coefficients, which allow to compute gamma up to the precision of either Float or Double correspondingly. These lists of coefficients were generated by LanczosConstants module.

If one will be able to run LanczosConstants for different floating types, he will be able to generate relevant coefficients and derive implementations of gamma for other types.

As you noted above, we do not need this at the moment for arithmoi. And I have not dig much into gamma. Just thinking aloud.

[rockbmb]

@Bodigrim sorry, my question was more about: you said something about moving some code to matrix, but also release it as a user-facing module of gamma. I thought one would do either one or the other.

[Bodigrim]

I meant "migrate LanczosConstants [from an unreleased library] to matrix". Currently LanczosConstants relies on an unreleased library for matrix computations, which has prevented the original author to upload it on Hackage as a part of gamma.

[rockbmb]

Oh, I understand - you meant to migrate the matrix dependency the module relies on, not the module itself.

[rockbmb]

By the way, I'm awaiting a reply on gamma, but so for now I'll focus on closing this PR.

It'll be difficult to make the Hurwitz zeta function numerically stable on real inputs, that can be a patch for a different occasion. For now I'll keep the first argument integral. I'll push what I have now regardless for future reference, same as what I did with the Eisenstein integer gcd algorithm, and then revert it.

[rockbmb]

@Bodigrim I've been reading re-reading the section on bounding the error for the Hurwitz zeta function, Arb's documentation, and I think that at this point it is not worth it to implement error bounding any more complex than what is already there.

An idea I've been considering lately is to create implement ball arithmetic types in GHC (in the Prelude, alongside Double and Float, perhaps?), and then write a small set of C FFI bindings to the Arb library (or any other that offers the same functionality, although Arb seems very good and, for now, maintained) only for the functions arithmoi currently uses. If there's a function truncate :: Floating a => Arb_t -> a, where Arb_t is a Haskell wrapper over arb_t from the Arb library, then the end-user can choose what type to run the computations with.
What do you think about this idea?

In any case, this PR is good to go on my end.

[rockbmb]

@Bodigrim by the way, when you're working on your PRs, how do you generate documentation locally to see how it's rendered in Hackage? cabal haddock --internal?

I just noticed that the documentation I wrote for the Eisenstein integers module looks pretty bad, I should check that sort of stuff before PRs get merged.

[Bodigrim]

I use stack haddock.

In any case, this PR is good to go on my end.

Good. I'll take a look during the next week.

[Bodigrim]

Please rebase atop latest master.

What I am worried about is that benchmarks of zeta/beta functions (kinda "sum first 20 values of Riemann zeta function) became significantly slower. This may be connected not only to possibly worse convergence properties, but also to the fact that zetasOdd recycle previous values to compute next ones, saving time. I believe it may be worth to leave implementation of Riemann zeta function as is, and add a test checking that it matches Hurwitz zeta function.

[rockbmb]

@Bodigrim that's probably because I implemented zetaHurwitz as a function and not as a recurrence, such as what zetas is now. I'll see if I can turn it into one, if not I'll implement your suggestion.

EDIT: If I recall correctly, it's because betas only needs zetaHurwitz at even values, while zetas only needs it for odd values. This means:

• half of the values in the hypothetical zetaHurwitz :: [a] will not be used (and so will be discarded) the way it's currently used, for both betas/zetas, or
• it can calculate every other value alone, based on a starting integer (2 for betas, 3 for zetas).

[rockbmb]

@Bodigrim I have returned to fix the issue with zetaHurwitz. The precision seems to have been kept, and now that there's more reusing of previous values as was done in zetasOdd/betasEven, performance should not be dissimilar (hopefully in the same order of magnitude).

I'll try benchmarking it now, but I'll be slightly skeptical about the numbers (running WSL).

[rockbmb]

From master:

Benchmark criterion: RUNNING...
benchmarked Zeta/riemann zeta
time                 75.59 μs   (75.16 μs .. 76.07 μs)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 75.10 μs   (74.79 μs .. 75.36 μs)
std dev              956.2 ns   (776.1 ns .. 1.210 μs)

benchmarked Zeta/dirichlet beta
time                 116.5 μs   (116.0 μs .. 117.2 μs)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 117.2 μs   (116.8 μs .. 117.6 μs)
std dev              1.299 μs   (1.049 μs .. 1.690 μs)

---

Current code (hurwitz as a recurrence):

benchmarked Zeta/riemann zeta
time                 145.2 μs   (125.7 μs .. 163.7 μs)
                     0.923 R²   (0.897 R² .. 0.952 R²)
mean                 148.3 μs   (142.0 μs .. 155.6 μs)
std dev              23.41 μs   (19.90 μs .. 28.43 μs)
variance introduced by outliers: 80% (severely inflated)

benchmarked Zeta/dirichlet beta
time                 290.6 μs   (259.6 μs .. 320.9 μs)
                     0.929 R²   (0.891 R² .. 0.960 R²)
mean                 258.4 μs   (247.4 μs .. 272.0 μs)
std dev              39.40 μs   (34.22 μs .. 45.39 μs)
variance introduced by outliers: 79% (severely inflated)

---

Previous code (hurwitz as a function):

benchmarked Zeta/riemann zeta
time                 663.2 μs   (654.8 μs .. 674.8 μs)
                     0.998 R²   (0.997 R² .. 0.999 R²)
mean                 692.2 μs   (685.7 μs .. 706.4 μs)
std dev              33.06 μs   (18.84 μs .. 64.64 μs)
variance introduced by outliers: 27% (moderately inflated)

benchmarked Zeta/dirichlet beta
time                 1.287 ms   (1.282 ms .. 1.292 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 1.286 ms   (1.283 ms .. 1.290 ms)
std dev              12.37 μs   (9.008 μs .. 19.14 μs)

[rockbmb]

Went from 10x slower to 2x, I guess that's barely acceptable. There might be room for improvement, @Bodigrim can you review the PR once more?

[Bodigrim]

It is certainly acceptable, good job! Overall it looks good to me, but I am wondering: do you know what exactly improved performance 5x? My expectation was that Hurwitz zeta functions on consequent integers do not share much, so I am quite surprised.

[rockbmb]

@Bodigrim Here's a short explanation:

• From the 4.8.5 formula, we see that the only variables that vary as a result of the arguments passed to the Hurwitz zeta function are s, a, N, M. Since N ~ M and the algorithm I wrote is quite simple, using no clever ideas to find an appropriate N and just assuming N = M, what happens is that s, a, N are what we care about.
• For a particular kind of function, a will not change: if we're calculating a sequence of values of the Riemann zeta function at odd integers, a will be 1. If we're calculating a sequence of values of the Dirichlet beta function at even integers, a will be both 0.25 and 0.75 (this is why the benchmarks for Dirichlet beta are 2x as slow as Riemann's - it has to calculate two recurrences, one for a = 0.25 and one for a = 0.75).
• For the same desired precision, N will not change. If we pass 1e-15 to zetas or betas, then N's value will be the same for s=0, 1, 2, 3 ...

This means the only variable that iterates is s, giving the opportunity for a lot of sharing of previous values. For example, the fraction involving even indices of Bernoulli's recurrence divided by factorials of even indices was being recalculated every time zetaHurwitz was called. The same goes for rising factorials, and every term involving (a + N). Every such term would be calculated several times for each call of zetaHurwitz, but now it only happens once in the body of the zetaHurwitz recurrence.

[rockbmb]

I should also add that although there is still a performance penalty, there is a precision increase: see 4804b57. Doing this with the current code causes the code to fail, which means the current implementation offers higher precision for higher indices of the zetas/betas recurrences.

It might get even better after dmcclean/exact-pi#8 gets merged, although I don't think think it's going to happen for a while.

[rockbmb]

Remembered to fix Haddock, unlike with the EisensteinIntegers module that I only fixed after the fact.
@Bodigrim PR is good to go (unless you want to review it further).

[Bodigrim]

Thanks for awesome work!

[rockbmb]

@Bodigrim No problem. I should also thank you for taking so much of your time to review my PRs, give advice and merge code 👍 .