Zeta Functions & Combinatorics Example

Demonstrates riemann_zeta(), hurwitz_zeta(), expint(), and comb().

Example Code

"""Example: IMSL Special Functions — Riemann/Hurwitz Zeta, expint, and comb."""
from __future__ import annotations

import math
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from specialfunctions import riemann_zeta, hurwitz_zeta, comb, expint

# ---------------------------------------------------------------------------
# Known-value checks
# ---------------------------------------------------------------------------
print(f"riemann_zeta(2) = {riemann_zeta(2):.10f}  (expected: {math.pi**2/6:.10f}  = pi^2/6, Basel problem)")
print(f"riemann_zeta(4) = {riemann_zeta(4):.10f}  (expected: {math.pi**4/90:.10f}  = pi^4/90)")
print(f"hurwitz_zeta(2, 1) = {hurwitz_zeta(2, 1):.10f}  (expected: {math.pi**2/6:.10f}  = zeta(2))")
print(f"comb(10, 5)  = {comb(10, 5)}  (expected: 252)")
print(f"comb(6, 0)   = {comb(6, 0)}  (expected: 1)")
print(f"comb(6, 6)   = {comb(6, 6)}  (expected: 1)")
print(f"expint(1, 1) = {expint(1, 1):.10f}")
print(f"expint(1, 0) = {expint(1, 0):.6f}  (E1(0) = +inf)")

# ---------------------------------------------------------------------------
# Plotting
# ---------------------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle("Riemann ζ, Hurwitz ζ, Exponential Integral Eₙ, and Binomial Coefficients", fontsize=13)

# --- Panel 1: Riemann zeta ζ(s) for s in (1, 10] --------------------------
ax = axes[0, 0]
s_vals = np.linspace(1.05, 10, 500)
ax.plot(s_vals, riemann_zeta(s_vals), "tab:blue", lw=1.8, label="ζ(s)")
ax.axhline(1.0, color="k", lw=0.7, ls=":", label="ζ→1 as s→∞")
ax.axhline(math.pi**2 / 6, color="tab:orange", lw=1.2, ls="--",
           label=f"ζ(2)=π²/6≈{math.pi**2/6:.4f}")
ax.set_title("Riemann Zeta ζ(s)  (real axis, s > 1)")
ax.set_xlabel("s")
ax.set_ylabel("ζ(s)")
ax.set_ylim(0.9, 6)
ax.legend()
ax.grid(True, alpha=0.3)

# --- Panel 2: Hurwitz zeta ζ(2, a) vs a in (0, 3] -------------------------
ax = axes[0, 1]
a_vals = np.linspace(0.1, 4.0, 400)
ax.plot(a_vals, hurwitz_zeta(2, a_vals), "tab:green", lw=1.8, label="ζ(2, a)")
ax.axvline(1.0, color="k", lw=0.7, ls=":", label="a=1  (Riemann ζ(2))")
ax.set_title("Hurwitz Zeta ζ(2, a) vs a")
ax.set_xlabel("a")
ax.set_ylabel("ζ(2, a)")
ax.set_ylim(0, 15)
ax.legend()
ax.grid(True, alpha=0.3)

# --- Panel 3: Exponential integral E_n(x) for n = 1, 2, 3 -----------------
ax = axes[1, 0]
x_ei = np.linspace(0.05, 4, 400)
colours = ["tab:blue", "tab:orange", "tab:green"]
for n, colour in zip([1, 2, 3], colours):
    ax.plot(x_ei, expint(n, x_ei), color=colour, lw=1.8, label=f"E{n}(x)")
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.set_title("Generalised Exponential Integral Eₙ(x)")
ax.set_xlabel("x")
ax.set_ylabel("Eₙ(x)")
ax.set_ylim(0, 3)
ax.legend()
ax.grid(True, alpha=0.3)

# --- Panel 4: Pascal's triangle using comb(n, k) ---------------------------
ax = axes[1, 1]
n_rows = 9
for row in range(n_rows):
    for k in range(row + 1):
        val = comb(row, k)
        x_pos = k - row / 2
        y_pos = -row
        ax.text(x_pos, y_pos, str(val), ha="center", va="center",
                fontsize=8, bbox=dict(boxstyle="round,pad=0.2", fc="lightyellow", ec="gray", lw=0.5))
ax.set_xlim(-n_rows / 2 - 0.5, n_rows / 2 + 0.5)
ax.set_ylim(-n_rows + 0.3, 1)
ax.set_title("Pascal's Triangle — comb(n, k)")
ax.axis("off")

plt.tight_layout()
plt.savefig("test_output/example_imsl_zeta_combinatorics.svg", bbox_inches="tight")
print("Saved: test_output/example_imsl_zeta_combinatorics.svg")

Plot Output

Riemann ζ(s), Hurwitz ζ(2,a), generalised exponential integral Eₙ(x), and Pascal’s triangle generated via comb(n,k).

Console Output

riemann_zeta(2) = 1.6449340668  (expected: 1.6449340668  = pi^2/6, Basel problem)
riemann_zeta(4) = 1.0823232337  (expected: 1.0823232337  = pi^4/90)
hurwitz_zeta(2, 1) = 1.6449340668  (expected: 1.6449340668  = zeta(2))
comb(10, 5)  = 252  (expected: 252)
comb(6, 0)   = 1  (expected: 1)
comb(6, 6)   = 1  (expected: 1)
expint(1, 1) = 0.2193839344
expint(1, 0) = inf  (E1(0) = +inf)
Saved: test_output/example_imsl_zeta_combinatorics.svg