Am Donnerstag, den 21.01.2010, 13:19 -0500 schrieb Raymond Toy:
> 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)));
Hello Ray,
yes, I think we should implement this algorithm. I think it is a
variation of an algorithm presented first in a paper of P. Borwein "An
Efficient Algorithm for the Riemann Zeta Function". In this paper it is
the third algorithm. It is a easier to implement then the second
algorithm presented in the paper of Borwein.
The algorithm will work for complex values too.
Do you already have an implementation using the bigfloat package?
Dieter Kaiser