stdlib
stdlib copied to clipboard
stdlib-js
[RFC]: Why not use Riemann-Siegel Formulation to very quickly evaluate zeta(s) for real part 1/2?
Description
Since the riemann zeta function is very interesting in the critical strip of Re(s)=1/2 I think it would be very useful to implement a very fast method to evaluate zeta on that strip. The Riemann-Siegel Formula is one that does exactly that. The formula is not simple but I have worked it out and implemented it in js before. It is as follows
$\zeta(\frac{1}{2} + it) = Z(t)e^{-i\vartheta(t)}$ where $\vartheta(t)=\Im[\log \Pi(\frac{it}{2}-\frac{3}{4})] - \frac{t}{2}\log 2\pi$ I have derived a very nice asymptotic series for this function using the Stirling Series of the Pi (factorial) function $\vartheta(t)=\frac{1}{2}t\bigg(\log\bigg(\frac{t}{2\pi}\bigg)-1\bigg) - \frac{\pi}{8}+\displaystyle \sum_{k=1}^\infty \frac{(1-2^{1-2k})|B_{2k}|}{4k(2k-1)t^{2k-1}}+\frac{1}{2}\arctan(e^{-\pi t})$ For the first few terms we have $\vartheta(t)=\frac{t}{2}\log(\frac{t}{2\pi}) - \frac{t}{2} - \frac{\pi}{8} + \frac{1}{48t} + \frac{7}{5760t^3} + \frac{31}{80640t^5} ... $ Next we tackle the main function: $Z(t)$ It is defined to be $Z(t)=\displaystyle 2\sum_{n=1}^N n^{-1/2} \cos(\vartheta(t)-t\log n) + R$ Where $N = \lfloor (\frac{t}{2\pi})^{1/2}\rfloor$ It is fairly simple but the main difficulty comes from the calculation of the remainder term R $R = (-1)^{N-1} (\frac{t}{2\pi})^{-1/4}[C_0 + C_1 (\frac{t}{2\pi})^{-1/2} + C_2 (\frac{t}{2\pi})^{-2/2} + C_3 (\frac{t}{2\pi})^{-3/2} + C_4 (\frac{t}{2\pi})^{-4/2} ...]$ The coefficients $C_n$ have horrificly large expansions. They are defined in terms of the $\Psi(x)$ function: $\Psi(x)=\frac{\cos( 2\pi(x^2-x-1/16))}{\cos(2\pi x)}$ and $p=(\frac{t}{2\pi})^{1/2} - N$ $C_0 = \Psi(p)$ $C_1 = -\frac{1}{96 \pi^2} \Psi^{(3)}(p)$ $C_2 = \frac{1}{18432 \pi^4} \Psi^{(6)}(p) + \frac{1}{64 \pi^2} \Psi^{(2)}(p)$ $C_3 = -\frac{1}{5308416 \pi^6} \Psi^{(9)}(p) - \frac{1}{3804 \pi^4} \Psi^{(5)}(p) - \frac{1}{64 \pi^2} \Psi^{(1)}(p)$ $C_4 = \frac{1}{2038431744 \pi^8} \Psi^{(12)}(p) + \frac{11}{5898240 \pi^6} \Psi^{(8)}(p) + \frac{19}{24576 \pi^4} \Psi^{(4)}(p)+\frac{1}{128 \pi^2} \Psi(p)$ As you can see the high derivatives of $\Psi(x)$ will yield massive expressions. Therefore it is important to use some sort of approximation for them. In my implementation I used a Chebyshev Approximation and used that as a polynomial array to take multiple derivative of the function. This is the hardest part of the formulation. The calculation of those random coefficients inside the $C_n$ terms is roughly equivalent to hell to compute and required a symbolic math system. I used sympy to calculate mine but its not needed to have more than 4 terms. If you think this would be a useful addition I can help you guys implement it. Thanks
Related Issues
No response
Questions
No.
Other
No.
Checklist
- [X] I have read and understood the Code of Conduct.
- [X] Searched for existing issues and pull requests.
- [X] The issue name begins with
RFC:.
:tada: Welcome! :tada:
And thank you for opening your first issue! We will get back to you shortly. :runner: :dash:
@Bobingstern Thanks for filing this issue. A couple questions:
- Are there precedents for the proposed function in other numerical libraries/languages? E.g., NumPy, SciPy, MATLAB, Julia, R, etc. If so, what is the API (signature), and what are the reference implementations?
- How does the proposed function relate to the existing stdlib Riemann Zeta function?
Mathematica, and MATLAB compute the zeta function over the entire complex plane (a+ib).
https://www.mathworks.com/help/symbolic/sym.zeta.html
https://reference.wolfram.com/language/ref/Zeta.html
The formula is connected to the existing zeta implementation by an extension. The library computes the zeta function for real values and this formula computes it for complex values with real part 1/2 (the critical line).
I am interested to work on this issue under Quine. So please assign me this task to me.
They only thing that is difficult about this formulation is the computation of $\Psi(x)$ and it's derivatives. If you have a clever way to go about it that would be fantastic!