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

enter image description here

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


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.


  • 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


2, 0


1.6449340668482266, 0


1, 1


0.5821580597520037, -0.9268485643308071


-1, 0


-0.08333333333333559, 0

  • 1
    \$\begingroup\$ What's the required output precision? I'm not sure I understand Your code must contain value that gives theoretically arbitrary precision when made arbitrary large/small. Do you mean like a loop maximum value than when increased without bound gives increased precision? Can that value be hard-coded? \$\endgroup\$ – Luis Mendo Mar 10 '16 at 16:32
  • \$\begingroup\$ @DonMuesli This means the precision depends on a parameter, say N, which you may give any value you like, but for any given precision, you can make N small or large enough to achieve that precision. The word theoretically is there because you must not worry about the limited precision of the machine or language. \$\endgroup\$ – Jens Renders Mar 10 '16 at 16:35
  • \$\begingroup\$ To further clarify N: is it sufficient that for any bound eps and input x there exists an N which calculates zeta(x) to within eps; or must there exist an N which depends only on eps and guarantees that for any x (or perhaps for any x more than a given function of eps from the pole) it achieves the bound; or may N depend on x, but answers should explain how to calculate N given x and eps? (My analytical number theory isn't up to much, but I suspect that options 2 and 3 are going to be beyond all but one or two regular posters). \$\endgroup\$ – Peter Taylor Mar 10 '16 at 18:59
  • \$\begingroup\$ @PeterTaylor N large enough: For any x and for any eps there must exist a P such that for all N>P the output is closer than eps to the exact value. Is this clear? Do I need to clarify it for the case with N small enough? \$\endgroup\$ – Jens Renders Mar 10 '16 at 19:27
  • \$\begingroup\$ No, that's clear enough. \$\endgroup\$ – Peter Taylor Mar 10 '16 at 20:40

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):
  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)
  if a*a+b*b<E:break
  • \$\begingroup\$ z(2,0) gives an incorrect value, should be pi^2/6. \$\endgroup\$ – GuillaumeDufay Sep 22 '18 at 2:13

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):
 while a*a+b*b>E:
  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)

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)    

     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])
  -- (m+iq)(r1+ir2)=(m*r1-q*r2)+i(m*r2+q*r1)

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)  []

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.