Alex Herbert created NUMBERS-216:
------------------------------------
Summary: Add a zeta function
Key: NUMBERS-216
URL: https://issues.apache.org/jira/browse/NUMBERS-216
Project: Commons Numbers
Issue Type: New Feature
Components: gamma
Reporter: Alex Herbert
Add a zeta function to evaluate the Hurwitz zeta function:
{noformat}
oo 1
zeta(s, a) = sum ------
k=0 s
(k+a)
{noformat}
The series is formally defined for complex s with real(s) > 1 and real a != 0,
-1, -2, ... It can be extended to any s != 1 using analytic continuation.
When a=1 this is the Riemann zeta function.
The Hurwitz zeta can be evaluated for real s != 1 and a > 0 efficiently using
the Euler-Maclaurin formula (see [Huwitz zeta function in DLMF 25.17 Eq
7|https://dlmf.nist.gov/25.11#E7]). The final integral can be dropped as a
residual error term. The formula uses a summation of n terms and is valid when
s > -2n and a > 0.
An implementation was added to Commons Statistics to support the Zipf and zeta
distributions. The implementation supports s > 1 and a > 0.
h2. Increasing support
Supporting large negative s requires long computation times as the required
number of terms n increases. Each term requires a Bernoulli number which can be
precomputed for small n but are expensive to compute dynamically.
Supporting negative a requires computing the summation of the formal series of
(k+a)^-s until k+a is above 0. There is no general reflection formula. This
will have long run times for negative a.
I suggest an initial implementation supporting only s > 1 where the series is
absolutely convergent. Negative a can be supported with the caveat of reduced
accuracy and long run times.
When a = 1 the function is the Riemann zeta function. There exists a reflection
formula for negative s allowing computation for any s != 1.
I have tested a double-precision implementation from [Boost C++
zeta|https://www.boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/zetas/zeta.html]
converted to Java using the same range supported by the current hurwitz zeta
implementation. Relative error on 5000 values with s in [1, 32):
||Function||Max error||RMS error||
|zeta(s, 1)|1.664744282142072|0.48639150510073753|
|zeta(s)|1.448659608140353|0.3295234039818026|
The Boost function has a small accuracy improvement and will be an order of
magnitude faster as it uses a polynomial approximation and avoids calls to
Math.pow. The Boost license is permissive and other Boost function
implementations already exist in Commons Numbers.
--
This message was sent by Atlassian Jira
(v8.20.10#820010)
