Bug 1995531
<https://sourceforge.net/tracker/?func=detail&aid=1995531&group_id=4933&atid=104933>;
mentions problems with zeta near 0.
I can confirm these, except zeta(0) no longer gives a division by zero
error. It correctly returns -1/2.
Here is an alternative implementation of bfzeta (in
share/numeric/bffac.mac). I don't know what algorithm is used there,
but, if nothing else, this new implementation is (gasp!) documented.
With this version, we have:
(%i132) rltzeta(1b-7);
(%o132) - 5.0000009189386335225b-1
(%i133) rltzeta(-1b-20);
(%o133) - 5.1145234580810622638b-1
(%i134) rltzeta(10b-13);
(%o134) - 5.0000000000091893854b-1
I haven't independently verified these results, but they make some
sense. The implementation needs some work because it only works with
bfloats now.
Ray
/*
*
* See http://numbers.computation.free.fr/Constants/constants.html
* and, in particular,
http://numbers.computation.free.fr/Constants/Miscellaneous/zetaevaluations.pdf.
* We use the algorithm from Proposition 2:
*
* zeta(s) = 1/(1-2^(1-s)) *
* (sum((-1)^(k-1)/k^s,k,1,n) +
1/2^n*sum((-1)^(k-1)*e(k-n)/k^s,k,n+1,2*n))
* + g(n,s)
*
* where e(k) = sum(binomial(n,j), j, k, n) and
* |g(n,s)| <= 1/8^n * 4^abs(sigma)/abs(1-2^(1-s))/abs(gamma(s))
*
* with s = sigma + %i*t
*
* This technique converge for sigma > -(n-1)
*
* Figure out how many terms we need from |g(n,s)|. The answer is
*
* n = log(f/eps)/log(8), where f =
4^abs(sigma)/abs(1-2^(1-s))/abs(gamma(s))
*/
rltzetahelp(n,k):=sum(bfloat(binomial(n,j)),j,k,n);
rltzeta(s):=
if s = 0 then -1/2
elseif s = 0.0 then -0.5
elseif s = 0b0 then -0.5b0
elseif is(s < 0) then
bfloat(2^s*%pi^(s-1)*sin(%pi/2*s)*gamma(1-s)) * rltzeta(1-s)
else
/*
* What should we do when s is a positive even integer for each an
exact
* value is known?
*/
block([n :
max(1,ceiling(log(4^abs(realpart(s))/abs(1-2^(1-s))/abs(gamma(s))/10^(-fpprec))/log(8)))],
print("n = ", n),
1/(1-2^(1-s)) *
(sum((-1)^(k-1)/bfloat(k^s),k,1,n) +
1/2^n*sum((-1)^(k-1)*rltzetahelp(n,k-n)/k^s,k,n+1,2*n)));