Evaluate the Riemann Zeta Function at a Complex Number

Question

Introduction

I found this question that was closed because it was unclear, yet it was a nice idea. I'll do my best to make this into a clear challenge.

The Riemann Zeta function is a special function that is defined as the analytical continuation of

to the complex plane. There are many equivalent formulas for it which makes it interesting for code golf.

Challenge

Write a program that takes 2 floats as input (the real and imaginary part of a complex number) and evaluates the Riemann Zeta function at that point.

Rules

• Input and output via console OR function input and return value
• Built in complex numbers are not allowed, use floats (number, double, ...)
• No mathematical functions except + - * / pow log and real valued trig functions (if you want to integrate, use the gamma function, ... you must include this functions definition in the code)
• Input: 2 floats
• Output: 2 floats
• Your code must contain value that gives theoretically arbitrary precision when made arbitrary large/small
• The behaviour at input 1 is not important (this is the only pole of this function)

Shortest code in bytes wins!

Example Input and Output

Input:

2, 0

Output:

1.6449340668482266, 0

Input:

1, 1

Output:

0.5821580597520037, -0.9268485643308071

Input:

-1, 0

Output:

-0.08333333333333559, 0

Answer 1

Python - 385

This is a straightforward implementation of Equation 21 from http://mathworld.wolfram.com/RiemannZetaFunction.html This uses Python's convention for optional arguments; if you want to specify a precision, you can pass a third argument to the function, otherwise it uses 1e-24 by default.

import numpy as N
def z(r,i,E=1e-24):
 R=0;I=0;n=0;
 while(True):
  a=0;b=0;m=2**(-n-1)
  for k in range(0,n+1):
   M=(-1)**k*N.product([x/(x-(n-k))for x in range(n-k+1,n+1)]);A=(k+1)**-r;t=-i*N.log(k+1);a+=M*A*N.cos(t);b+=M*A*N.sin(t)
  a*=m;b*=m;R+=a;I+=b;n+=1
  if a*a+b*b<E:break
 A=2**(1-r);t=-i*N.log(2);a=1-A*N.cos(t);b=-A*N.sin(t);d=a*a+b*b;a=a/d;b=-b/d
 print(R*a-I*b,R*b+I*a)

Answer 2

Python 3, 303 297 bytes

This answer is based on R.T.'s Python answer with several modifications:

• First, Binomial(n, k) is defined as p = p * (n-k) / (k+1) which changes Binomial(n,k) to Binomial(n,k+1) with every pass of the for loop.
• Second, (-1)**k * Binomial(n,k) became p = p * (k-n) / (k+1) which flips the sign at every step of the for loop.
• Third, the while loop has been changed to immediately check if a*a + b*b < E.
• Fourth, the bitwise not operator ~ is used in several places where they would aid in golfing, using identities such as -n-1 == ~n, n+1 == -~n, and n-1 == ~-n.

Several other small modifications were made for better golfing, such as putting the for loop on one line and the call to print on one line with the code before it.

Golfing suggestions welcome. Try it online!

Edit: -6 bytes from a number of small changes.

import math as N
def z(r,i,E=1e-40):
 R=I=n=0;a=b=1
 while a*a+b*b>E:
  a=b=0;p=1;m=2**~n
  for k in range(1,n+2):M=p/k**r;p*=(k-1-n)/k;t=-i*N.log(k);a+=M*N.cos(t);b+=M*N.sin(t)
  a*=m;b*=m;R+=a;I+=b;n+=1
 A=2**-~-r;t=-i*N.log(2);x=1-A*N.cos(t);y=A*N.sin(t);d=x*x+y*y;return(R*x-I*y)/d,(R*y+I*x)/d

Answer 3

Axiom, 413 315 292 bytes

p(n,a,b)==(x:=log(n^b);y:=n^a;[y*cos(x),y*sin(x)]);z(a,b)==(r:=[0.,0.];e:=10^-digits();t:=p(2,1-a,-b);y:=(1-t.1)^2+t.2^2;y=0=>[];m:=(1-t.1)/y;q:=t.2/y;n:=0;repeat(w:=2^(-n-1);abs(w)<e=>break;r:=r+w*reduce(+,[(-1)^k*binomial(n,k)*p(k+1,-a,-b) for k in 0..n]);n:=n+1);[r.1*m-q*r.2,m*r.2+r.1*q])

This would implement too the equation 21 from http://mathworld.wolfram.com/RiemannZetaFunction.html The above should be the one iterpreted Axiom function z(a,b) here 16x slower than this below function Zeta(a,b)[that should be the one compiled] all ungolfed and commented [1 second for Zeta() against 16 seconds for z() for one value of 20 digits afther the float point]. For the digit question, one would choose the precision by calling digits(); function, for example digits(10);z(1,1) should print 10 digits after the point, but digits(50);z(1,1) should print 50 digits after the point.

-- elevImm(n,a,b)=n^(a+i*b)=r+i*v=[r,v]
elevImm(n:INT,a:Float,b:Float):Vector Float==(x:=log(n^b);y:=n^a;[y*cos(x),y*sin(x)]::Vector Float);

--                      +oo               n
--                      ---              ---
--             1        \       1        \            n 
--zeta(s)= ---------- * /     ------  *  /    (-1)^k(   )(k+1)^(-s)
--          1-2^(1-s)   ---n  2^(n+1)    ---k         k  
--                       0                0


Zeta(a:Float,b:Float):List Float==
  r:Vector Float:=[0.,0.]; e:=10^-digits()

  -- 1/(1-2^(1-s))=1/(1-x-i*y)=(1-x+iy)/((1-x)^2+y^2)=(1-x)/((1-x)^2+y^2)+i*y/((1-x)^2+y^2)    

  t:=elevImm(2,1-a,-b);
  y:=(1-t.1)^2+t.2^2;
  y=0=>[] 
  m:=(1-t.1)/y; 
  q:=t.2/y
  n:=0
  repeat
     w:=2^(-n-1)
     abs(w)<e=>break  --- this always terminate because n increase
     r:=r+w*reduce(+,[(-1)^k*binomial(n,k)*elevImm(k+1,-a,-b) for k in 0..n])
     n:=n+1
  -- (m+iq)(r1+ir2)=(m*r1-q*r2)+i(m*r2+q*r1)
  [r.1*m-q*r.2,m*r.2+r.1*q]

this is one test for the z(a,b) function above:

(10) -> z(2,0)
   (10)  [1.6449340668 482264365,0.0]
                                              Type: List Expression Float
(11) -> z(1,1)
   (11)  [0.5821580597 520036482,- 0.9268485643 3080707654]
                                              Type: List Expression Float
(12) -> z(-1,0)
   (12)  [- 0.0833333333 3333333333 3,0.0]
                                              Type: List Expression Float
(13) -> z(1,0)
   (13)  []