[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](https://mathworld.wolfram.com/ChebyshevApproximationFormula.html) 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](https://github.com/stdlib-js/stdlib/blob/develop/CODE_OF_CONDUCT.md).
- [X] Searched for existing issues and pull requests.
- [X] The issue name begins with `RFC:`.