riemann.cpp

/**
   Compile:   g++ -O2 -g -o riemann riemann.cpp

   Run:       ./riemann 1.1 | tee riemann.dat
              (try values between 1 and 2)

   Gnuplot:   plot "riemann.dat" u 1:2 w lines
 */

#include <complex>
#include <cmath>
#include <cstdio>

typedef std::complex<double> cdouble;
static const cdouble j{0.0, 1.0};

cdouble zeta(double g, double t, double tol = 1e-3)
{
    cdouble z = 0.0;
    cdouble s = g + j * t;
    for (int n = 1; ; ++n) {
        cdouble term = 1.0 / std::pow(n, s);
        z += term;
        if (std::abs(term) < tol)
            break;
    }
    return z;
}

int main(int argc, char *argv[])
{
    if (argc != 2)
        return 1;
    double g = atof(argv[1]);

    const double x1 = 0.0;
    const double x2 = 50.0;
    const double xstep = 0.05;
    for (double x = x1; x <= x2; x += xstep) {
        cdouble z = zeta(g, x);
        printf("%f %f %f %f\n", x, std::abs(z), z.real(), z.imag());
    }

    return 0;
}