Programming details on the computation of the Riemann zeta function using Aribas

RamaWolf
Messages
95
Reaction score
2
(1) Let s be a complex number like s = a + b i, then we define \zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}

Our aim:

to compute ζ(\frac{1}{2}+14.1347 i) with the help of the programming language Aribas

(2) Web Links

Aribas: http://www.mathematik.uni-muenchen.de/~forster/sw/aribas.html

Dirichlet eta: http://en.wikipedia.org/wiki/Dirichlet_eta_function

Euler transformation of alterning series: http://en.wikipedia.org/wiki/Euler_transform

(3) Books: Knopp: Theorie und Anwendung der unendlichen Reihen; Henrici: Applied
and Computational Complex Analsis; Derbyshire: Prime Obsession

(4) Basics: The Riemann zeta in (1) is defined only for s.Re > 1 (i.e. the real part of s);
we use the 'alternating zeta' or Dirichlet eta defined as

\eta(s) = \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n^s}

defined for s.Re > 0

The Riemann zeta and the Dirichlet eta are related through:

ζ(s) := η(s) / (1-2^{1-s})

With the Euler transformation of the Dirichlet eta, the Riemann zeta looks like:

\zeta(s)=\frac{1}{1-2^{1-s}}<br /> \sum_{n=0}^\infty \frac {1}{2^{n+1}}<br /> \sum_{k=0}^n (-1)^k {n \choose k} (k+1)^{-s}


To be continued with the programming details
 
Physics news on Phys.org
Aribas is a computer program for number theoretical calculations.
It comes with:
- a built-in large integer arithmetic
- several modes for floating-point accuray; single/double/long (32/64/128 bit) and user-defined bit number
- a bunch of number theoretical functions
- a great variety of example program scripts
- the possibility to easy implement additional user-written functions

For the Riemann zeta calculations, we have to implement a complex arithmetic,
which is easily done by a handfull of functions; we need especially;

AddC(a,b) to add two complex variables, SubC(a,b), AbsC(a), which
can be programmed at a very basic programming level.

A little bit more is needed for a function to calculate the
exponentiation of a real base raised to a complex exponent.

As an example, here is the Aribas code for 'PowRC':
(a complex variable c is stored in an array[2] of reals,
the Re part in c[0], the Im part in c[1])

Code:
function PowRC(b: real; a:complex): complex;
var
   c: complex;
   x,y: real;
begin 
   x:=log(b);
   y:=x*a[1];
   x:=exp(a[0]*x);
   c[0]:=x*cos(y);
   c[1]:=x*sin(y); 
   return(c);
end.

To be continued with the further programming details
 
Code:
Deta(x,n)={
	my(b=2^(2*n-1), c=b, s=0);
	forstep(k=n-1,0,-1,
		s += c*(-1)^(k)/((k+1)^x);
		b *= ((2*k+1)*(k+1)) / (2*(n+k)*(n-k));
		c += b;
	);
	return(s/c);
}
Rzeta(s,n=100)=Deta(s,n)/(1-2^(1-s))
 
Aribas code for the Dirichlet eta function:

Code:
   i:=1; t:=(0.0,0.0);
   for n:=0 to 9999 do
      i:=i*2;      
      x:=(0.0,0.0);
      for k:=0 to n do
         y:=PowRC(1.0/(1.0+k),s);
         j:=nCr(n,k); 
         y:=j*y; 
         if odd(k) then
            x:=SubC(x,y);
         else
            x:=AddC(x,y);
         end;
      end;
      x:=x/i;
      t:=AddC(t,x);
      if AbsC(x) < eps then
         [B]break[/B];
      end;
   end;
   return(t);


The for loop will terminate after 9999 iterrations, but in practice, the break
will be effective miuch earlier

{n \choose r} the binomial coefficient is computed with the function 'nCr' ("from n choose r")

Numerical results for s = \frac{1}{2}+14.1347 I

eps = 1.00000000000000000e-13

eta2: -0.000002874863 - 0.000047246855 I

factor: 0.411257843979 + 0.091385435419 I

RiemannZeta(s): 0.000003135364 - 0.000019693360 I
 
realprecision = 28 significant digits
eps=1e-13
zeta embeded function in PARI/GP
Code:
s=1/2.+I*14.1347
abs(YourZeta(s)-zeta(s))==2.4718664606232304313 E-14  60 iterations 219 mSec
abs(MyZeta(s)  -zeta(s))==1.6630931417336143861 E-15  30 iterations  15 mSec  
abs(zeta(s)    -zeta(s))==0.E-33                      ??              0 mSec
 
Last edited:
Back
Top