Riemann Zeta Function (C++)

The Riemann Zeta Function, denoted as ζ(s), is a mathematical function that takes a complex number s as input and sums an infinite series. For a basic understanding, when s is a real number greater than 1, it can be written as:

ζ(s) = 1⁻ˢ + 2⁻ˢ + 3⁻ˢ + 4⁻ˢ + ...

‍

In simpler terms, it adds up the reciprocals of all positive integers raised to the power s. For example:

• If s = 2, then ζ(2) = 1/1² + 1/2² + 1/3² + ... = 1 + 1/4 + 1/9 + ... ≈ 1.64493 (this is π²/6).

‍

However, in the case of Riemann, we are interested in only complex numbers. It is perfectly valid for real numbers but it only reveals its power in the complex plane. In fact, it is most interesting where the real part of the complex number is greater than zero but less than one. Known as the "critical strip"

‍

The code below focuses on an even more specific case where the function has real part 1/2 that will perhaps one day become obvious as to why. If you do run this code try 10,000 terms and choose imaginary part 14.134725 . It should get you close to zero.

‍

#include <iostream>
#include <complex>
#include <cmath>

using namespace std;

// Dirichlet series for Zeta (Re(s) > 1)
complex<double> zetaDirichlet(complex<double> s, int terms) {
    // Note: Simple sum for Re(s) > 1
    complex<double> sum(0.0, 0.0);
    for (int n = 1; n <= terms; n++) {
        sum += pow(1.0 / static_cast<double>(n), s);
    }
    return sum;
}

// Alternating series for Zeta in critical strip (Re(s) > 0)
complex<double> zetaAlternating(complex<double> s, int terms) {
    // Note: Uses ζ(s) = (1 - 2^(1-s))^(-1) * Σ ((-1)^(n-1) / n^s)
    // Converges for Re(s) > 0, including critical line Re(s) = 1/2
    complex<double> sum(0.0, 0.0);
    for (int n = 1; n <= terms; n++) {
        double sign = (n % 2 == 0) ? -1.0 : 1.0; // Alternates (-1)^(n-1)
        sum += sign * pow(1.0 / static_cast<double>(n), s);
    }
    complex<double> factor = 1.0 / (1.0 - pow(2.0, 1.0 - s));
    return factor * sum;
}

// Riemann Zeta Function
complex<double> riemannZeta(complex<double> s, int terms) {
    double re_s = real(s);
    if (re_s > 1.0) {
        // Note: Use Dirichlet series for Re(s) > 1
        return zetaDirichlet(s, terms);
    } else if (re_s > 0.0) {
        // Note: Use alternating series for 0 < Re(s) <= 1 (critical strip)
        return zetaAlternating(s, terms);
    } else {
        // Note: For Re(s) <= 0, we’d need functional equation, but not needed here
        cout << "Re(s) <= 0 not implemented in this demo." << endl;
        return complex<double>(0.0, 0.0);
    }
}

int main() {
    int terms;
    double t;

    cout << "Enter the number of terms for approximation: ";
    cin >> terms;
    cout << "Enter the imaginary part t (for s = 1/2 + it): ";
    cin >> t;

    // Note: s = 0.5 + it is on the critical line
    complex<double> s(0.5, t);
    complex<double> zeta = riemannZeta(s, terms);

    cout << "s = " << s << " (Re(s) = " << real(s) << ", Im(s) = " << imag(s) << ")" << endl;
    cout << "Zeta(s) approximation with " << terms << " terms: " << zeta << endl;
    cout << "Magnitude: " << abs(zeta) << endl;

    // Note: Check if close to a zero
    if (abs(zeta) < 0.001) {
        cout << "This is very close to a non-trivial zero!" << endl;
    }

    return 0;
}