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:
##\textbf{Exercise 10}:## I came across the following solution online: Questions: 1. When the author states in "that ring (not sure if he is referring to ##R## or ##R/\mathfrak{p}##, but I am guessing the later) ##x_n x_{n+1}=0## for all odd $n$ and ##x_{n+1}## is invertible, so that ##x_n=0##" 2. How does ##x_nx_{n+1}=0## implies that ##x_{n+1}## is invertible and ##x_n=0##. I mean if the quotient ring ##R/\mathfrak{p}## is an integral domain, and ##x_{n+1}## is invertible then...
The following are taken from the two sources, 1) from this online page and the book An Introduction to Module Theory by: Ibrahim Assem, Flavio U. Coelho. In the Abelian Categories chapter in the module theory text on page 157, right after presenting IV.2.21 Definition, the authors states "Image and coimage may or may not exist, but if they do, then they are unique up to isomorphism (because so are kernels and cokernels). Also in the reference url page above, the authors present two...
When decomposing a representation ##\rho## of a finite group ##G## into irreducible representations, we can find the number of times the representation contains a particular irrep ##\rho_0## through the character inner product $$ \langle \chi, \chi_0\rangle = \frac{1}{|G|} \sum_{g\in G} \chi(g) \chi_0(g)^*$$ where ##\chi## and ##\chi_0## are the characters of ##\rho## and ##\rho_0##, respectively. Since all group elements in the same conjugacy class have the same characters, this may be...
Back
Top