Connecting the Riemann Zeta Function with the Dirichlet Eta and Beta Functions

and Why Even Numbers of Zeta function Have Modulus $> 0$ Off the Critical Line

Author: Dominika Pillerová

Link:

• Pillerova, Dominika (2026): "Connecting the Riemann Zeta Function with the Dirichlet Eta and Beta Functions" Link: https://qiskit.cz/zeta-eta-beta-en.html

This notebook develops a unified framework connecting the Riemann zeta, Dirichlet eta, and Dirichlet beta functions through a single finite sine-weighted sum.

Key contributions:

1. Even-number representation of η and ζ. A finite sum of the form $\sum (n^s-1)/n^s \cdot \sin(\pi(1-n)/2)$ acts as an exact sieve, selecting only even $n$ and producing the eta function structure. After normalisation by $(1-2^{1-s})^{-1}$ and the exact scaling coefficient $C(s) = (\pi/2)^{1-s}$, the sum equals $\zeta(s)$ exactly. This shows that the zeros of ζ on the critical line are determined exclusively by even numbers.

2. Odd-number representation and the beta function. Shifting the sine argument by 1 — replacing $\sin(\pi(1-n)/2)$ with $\sin(\pi(0-n)/2)$ — switches the sieve to odd $n$, yielding the Dirichlet beta function $\beta(s)$, whose special value $\beta(2) = G$ (Catalan's constant) stands in contrast to the $\pi$-based values of $\zeta$ at even arguments.

3. The δ-function and rotation between β and ζ. A generalised function $\delta(k, s, N)$ with continuous parameter $k \in [0,1]$ interpolates between $\beta$ ($k=0$) and $\zeta$ ($k=1$) via the rotation identity $\delta(k,s) = \beta(s)\cos(\pi k/2) + (\eta(s)/2^s)\sin(\pi k/2)$. At $s=2$ the pair $(G,\, \pi^2/48)$ forms a 2D vector rotated by integer-entry matrix $M$, and the identity $\mathrm{Li}_2(i) = -\pi^2/48 + iG$ is verified to 10 decimal places.

$$\boxed{\delta(k, s, N) = \sum_{n=1}^{N} 2\cdot\frac{n^s-1}{n^s} \cdot\sin!\left(\frac{\pi(k-n)}{2}\right) \cdot\left(\frac{\pi}{2}\right)^{!1-s}}$$

$$\boxed{ \begin{pmatrix} \delta(0,2) \\ \delta(1,2) \end{pmatrix} = \underbrace{\begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}}_{M} \begin{pmatrix} G \\ \pi^2/48 \end{pmatrix} = \begin{pmatrix} -\pi^2/48 \\ G \end{pmatrix} }$$

$$\boxed{\delta(k, s, \infty) = \beta(s)\cdot\cos\!\left(\frac{\pi k}{2}\right) + \frac{\eta(s)}{2^s}\cdot\sin\!\left(\frac{\pi k}{2}\right)}$$

4. Why zeros cannot exist off the critical line for even numbers. The unnormalised $\delta_\text{even}$ satisfies $\delta_\text{even}(s) = \pi^{s-1}[\eta(s) - 2^{s-1}]$. At any zero $s_0$ of $\zeta$, $\eta(s_0)=0$ but $|2^{s_0-1}| = 1/\sqrt{2} \neq 0$, so $|\delta_\text{even}(s_0)| > 0$ off the critical line. The modulus is minimised exactly at $\sigma = 1/2$ due to the weight symmetry $n^{-\sigma} \leftrightarrow n^{-(1-\sigma)}$ imposed by the functional equation — confirming that even-number zeros are structurally confined to the critical line.

1. Classical Definition

The Riemann zeta function is defined for real $s > 1$ by the absolutely convergent series:

$$\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}$$

This definition fails for $s \leq 1$ — the series diverges. The goal is to understand how to extend this function to the entire complex plane (with the exception of the pole at $s = 1$).

2. Motivation: Why Does the Complex Plane Matter?

The Riemann Hypothesis asserts that all non-trivial zeros of $\zeta(s)$ lie on the line $\text{Re}(s) = \frac{1}{2}$. This line lies outside the region of convergence of the series $\sum n^{-s}$.

To speak about zeros at all, we need analytic continuation — the unique holomorphic extension of $\zeta$ from the region $\text{Re}(s) > 1$ to $\mathbb{C} \setminus \{1\}$.

3. Reformulation: Alternative Series Representations

As a first step, we examine algebraic rewritings of the classical series. The expression

$$\frac{n^s - 1}{n^s} = 1 - \frac{1}{n^s}$$

gives the identity:

$$N - \sum_{n=1}^{N} \frac{n^s - 1}{n^s} = \sum_{n=1}^{N} \frac{1}{n^s}$$

which is precisely $\zeta_N(s)$. This is a trivial rewriting — we verify numerically that both forms yield identical values.

Note: The sum $\sum (1 - n^{-s})$ alone diverges for all $s > 0$ and does not constitute a regularization in itself. It acquires meaning only after subtracting the divergent part $N$.

4. Numerical Exploration on the Real Line

We first visualize the behavior of $\zeta(s)$ for real $s \in (0, 20]$:

• for large $s$: $\zeta(s) \to 1$ (the $n=1$ term dominates)
• for $s = 2$: $\zeta(2) = \pi^2/6 \approx 1.6449$ (Euler's result)
• for $s \to 1^+$: $\zeta(s) \to +\infty$ (the harmonic series)

This is the starting point before extending to $\mathbb{C}$.

5. Comparison of Formulations

The plot below shows:

• Classical zeta $\zeta(s) = \sum n^{-s}$ (converges for $s > 1$, approximated here with 1000 terms)
• Rewritten form $N - \sum (n^s-1)/n^s$ — numerically identical, verifying algebraic equivalence

The agreement confirms consistency of both representations before we proceed to analytic continuation.

import numpy as np
import matplotlib.pyplot as plt

def modified_zeta(s, n_terms=10000):
    total = 0
    for n in range(1, n_terms + 1):
        total += ((n ** s) - 1) / (n ** s) 
    total = n_terms - total
    return total

def modified_zeta2(s, n_terms=10000):
    total = 0
    for n in range(1, n_terms + 1):
        total += ((n ** s) - 1) / (n ** s) 
    #total = n_terms - total
    return total

def classical_zeta(s, n_terms=10000):
    total = 0
    for n in range(1, n_terms + 1):
        total += 1 / (n ** s)
    return total

# Create s values (exponent)
s_values = np.linspace(0.01, 20, 500)

# Calculate both functions
zeta_modified = [modified_zeta(s, n_terms=1000) for s in s_values]
zeta_modified2 = [modified_zeta2(s, n_terms=1000) for s in s_values]
zeta_classical = [classical_zeta(s, n_terms=1000) for s in s_values]

# Define pastel colors
pastel_coral = '#FF6B6B'
pastel_blue = '#6BB6FF'

# Create plot
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 10))

# Plot 1: Modified zeta function
ax1.plot(s_values, zeta_modified, color=pastel_coral, linestyle='--', linewidth=2.5, 
         label='ζ(s) = N - Σ(n^s-1)/n^s')
ax1.plot(s_values, zeta_modified2, color=pastel_coral, linestyle='-', linewidth=2.5, 
         label='ζ(s) = Σ(n^s-1)/n^s')
ax1.grid(True, alpha=0.3)
ax1.set_xlabel('s', fontsize=14)
ax1.set_ylabel('ζ(s)', fontsize=14)
ax1.set_title('Modified Riemann Zeta Function: ζ(s) = Σ(n^s-1)/n^s', 
              fontsize=16, fontweight='bold')
ax1.legend(fontsize=12)
ax1.axhline(y=0, color='black', linestyle='-', linewidth=0.5)
ax1.set_ylim(bottom=0)

# Plot 2: Comparison of both functions
ax2.plot(s_values, zeta_modified, color=pastel_coral, linestyle='--', linewidth=2.5, 
         label='Modified: Σ(n^s-1)/n^s', alpha=0.8)
ax2.plot(s_values, zeta_classical, color=pastel_blue, linestyle='-', linewidth=2.5, 
         label='Classical: Σ1/n^s', alpha=0.8)
ax2.grid(True, alpha=0.3)
ax2.set_xlabel('s', fontsize=14)
ax2.set_ylabel('ζ(s)', fontsize=14)
ax2.set_title('Comparison: Modified vs. Classical Riemann Zeta Function', 
              fontsize=16, fontweight='bold')
ax2.legend(fontsize=12)
ax2.axhline(y=0, color='black', linestyle='-', linewidth=0.5)
ax2.set_ylim(bottom=0)

plt.tight_layout()
plt.show()

# Print some values
print("Modified Riemann Zeta Function: ζ(s) = Σ(n^s-1)/n^s")
print("=" * 60)
print(f"{'s':<10} {'ζ_modified(s)':<20} {'ζ_classical(s)':<20}")
print("-" * 60)

for s in [0.05, 0.07, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0]:
    zm = modified_zeta(s, n_terms=10000)
    zc = classical_zeta(s, n_terms=10000)
    print(f"{s:<10.2f} {zm:<20.6f} {zc:<20.6f}")

print("\nObservations:")
print(f"- For s=2: ζ_modified ≈ {modified_zeta(2, n_terms=10000):.6f}")
print(f"- For s=2: ζ_classical = π²/6 ≈ {np.pi**2/6:.6f}")

Modified Riemann Zeta Function: ζ(s) = Σ(n^s-1)/n^s
============================================================
s          ζ_modified(s)        ζ_classical(s)      
------------------------------------------------------------
0.05       6641.423150          6641.423150         
0.07       5642.783763          5642.783763         
0.10       4423.009022          4423.009022         
0.20       1980.461814          1980.461814         
0.30       900.494623           900.494623          
0.40       417.525500           417.525500          
0.50       198.544645           198.544645          
0.60       97.576122            97.576122           
0.70       50.052177            50.052177           
0.80       27.110644            27.110644           
0.90       15.688876            15.688876           
1.00       9.787606             9.787606            
1.50       2.592376             2.592376            
2.00       1.644834             1.644834            
2.50       1.341487             1.341487            
3.00       1.202057             1.202057            
4.00       1.082323             1.082323            
5.00       1.036928             1.036928            
6.00       1.017343             1.017343            
7.00       1.008349             1.008349            
8.00       1.004077             1.004077            
9.00       1.002008             1.002008            
10.00      1.000995             1.000995            
11.00      1.000494             1.000494            
12.00      1.000246             1.000246            
13.00      1.000123             1.000123            
14.00      1.000061             1.000061            
15.00      1.000031             1.000031            
16.00      1.000015             1.000015            
17.00      1.000008             1.000008            
18.00      1.000004             1.000004            
19.00      1.000002             1.000002            
20.00      1.000001             1.000001            

Observations:
- For s=2: ζ_modified ≈ 1.644834
- For s=2: ζ_classical = π²/6 ≈ 1.644934

6. Analytic Continuation via the Functional Equation

The key tool for extending $\zeta$ beyond the region $\text{Re}(s) > 1$ is the Riemann functional equation:

$$\zeta(s) = 2^s \pi^{s-1} \sin\!\left(\frac{\pi s}{2}\right) \Gamma(1-s)\, \zeta(1-s)$$

This identity, valid for all $s \in \mathbb{C} \setminus \{1\}$, allows us to evaluate $\zeta(s)$ for $s < 1$ using values of $\zeta(1-s)$, where $1-s > 0$ — that is, within the region of convergence of the original series.

We implement it as functional_equation_zeta(s).

7. Trivial Zeros and Special Values

The functional equation reveals the structure of $\zeta$ on the negative real axis. The factor $\sin(\pi s/2)$ vanishes at $s = -2, -4, -6, \ldots$ — these are the trivial zeros.

Special values arising from analytic continuation:

$s$	$\zeta(s)$	Note
$-1$	$-\frac{1}{12}$	Formally "$1+2+3+\ldots$"
$0$	$-\frac{1}{2}$
$2$	$\frac{\pi^2}{6}$	Basel problem (Euler)
$-2,-4,-6,\ldots$	$0$	Trivial zeros

Values such as $\zeta(-1) = -1/12$ are not the result of summing a divergent series — they are values of the analytically continued function at points where the original series is not defined at all.

8. Comparison: Truncated Sum vs. Analytic Continuation

The plot below compares:

• Analytic continuation (blue) — exact values of $\zeta(s)$ via the functional equation
• Truncated sum (red, dashed) — $\sum_{n=1}^{N} n^{-s}$ for finite $N$

For $s > 1$, the truncated sum approaches the analytic continuation (the larger $N$, the better).
For $s \leq 1$, the truncated sum diverges — the red curve does not represent $\zeta(s)$ there, only an artefact of finite $N$. Analytic continuation is the only correct definition of $\zeta$ in this region.

9. Summary: Strategy for Extension to the Complex Plane

The classical series $\sum n^{-s}$ converges only for $\text{Re}(s) > 1$.
For the entire complex plane (except $s = 1$) we will use in subsequent cells:

1. Dirichlet eta function $\eta(s) = \sum (-1)^{n+1} n^{-s}$ — converges for $\text{Re}(s) > 0$
2. Functional equation — covers the remainder of the complex plane
3. Euler–Maclaurin formula — efficient numerical approximation

The next step is the transition to the complex argument $s = \sigma + it$.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import zeta, gamma
import warnings
warnings.filterwarnings('ignore')

def modified_zeta_convergent(s, n_terms=10000):
    if s <= 0:
        return np.nan
    
    total = 0
    for n in range(1, n_terms + 1):
        #total += abs(1 - (n ** s)) / (n ** s)
        total += ((n ** s) - 1) / (n ** s)
    total = n_terms - total 
    
    return total

def functional_equation_zeta(s):
    if s >= 1:
        return zeta(s, 1)
    
    # Functional equation
    try:
        result = (2**s * np.pi**(s-1) * 
                 np.sin(np.pi * s / 2) * 
                 gamma(1 - s) * 
                 zeta(1 - s, 1))
        return result
    except:
        return np.nan

def zeta_analytical(s):
    """
    Complete analytical continuation for all s ≠ 1
    """
    if abs(s - 1) < 0.01:  # Near the pole
        return np.nan
    
    if s > 1:
        return zeta(s, 1)
    else:
        return functional_equation_zeta(s)

# Define pastel colors
pastel_coral = '#FF6B6B'
pastel_blue = '#6BB6FF'

# Create plot for real s values
s_values = np.linspace(-10, 10, 1000)

# Our modification (only for s > 0)
modified_values = []
for s in s_values:
    if s > 0.1 and s < 0.99:  # Avoid s=1
        modified_values.append(modified_zeta_convergent(s, n_terms=1000))
    elif s >= 1.01 and s <= 10:
        modified_values.append(modified_zeta_convergent(s, n_terms=5000))
    else:
        modified_values.append(np.nan)

# Analytical continuation
analytical_values = [zeta_analytical(s) for s in s_values]

# Create plots
fig, axes = plt.subplots(3, 1, figsize=(16, 14))

# Plot 1: Overall view
ax1 = axes[0]
ax1.plot(s_values, analytical_values, color=pastel_blue, linestyle='-', linewidth=2.5, 
         label='Analytical continuation ζ(s)', alpha=0.8)
ax1.plot(s_values, modified_values, color=pastel_coral, linestyle='--', linewidth=2, 
         label='My modification (s > 0)', alpha=0.8)
ax1.axhline(y=0, color='black', linestyle='-', linewidth=0.5, alpha=0.5)
ax1.axvline(x=0, color='black', linestyle='-', linewidth=0.5, alpha=0.5)
ax1.axvline(x=1, color='gray', linestyle='--', linewidth=1, alpha=0.5, label='s=1 (pole)')
ax1.grid(True, alpha=0.3)
ax1.set_xlabel('s', fontsize=14)
ax1.set_ylabel('ζ(s)', fontsize=14)
ax1.set_title('Analytical Continuation of Riemann Zeta Function', fontsize=16, fontweight='bold')
ax1.legend(fontsize=12)
ax1.set_ylim(-5, 10)

# Mark trivial zeros
trivial_zeros = [-2, -4, -6, -8, -10]
for zero in trivial_zeros:
    if -10 <= zero <= 10:
        ax1.plot(zero, 0, 'go', markersize=8, label='Trivial zeros' if zero == -2 else '')

# Plot 2: Detail for s > 0
ax2 = axes[1]
s_positive = s_values[s_values > 0]
analytical_positive = [zeta_analytical(s) for s in s_positive]
modified_positive = []
for s in s_positive:
    if s < 0.99 or s > 1.01:
        if s <= 5:
            modified_positive.append(modified_zeta_convergent(s, n_terms=5000))
        else:
            modified_positive.append(modified_zeta_convergent(s, n_terms=1000))
    else:
        modified_positive.append(np.nan)

ax2.plot(s_positive, analytical_positive, color=pastel_blue, linestyle='-', linewidth=2.5, 
         label='Analytical continuation', alpha=0.8)
ax2.plot(s_positive, modified_positive, color=pastel_coral, linestyle='--', linewidth=2, 
         label='My modification', alpha=0.8)
ax2.axhline(y=0, color='black', linestyle='-', linewidth=0.5, alpha=0.5)
ax2.axvline(x=1, color='gray', linestyle='--', linewidth=1, alpha=0.5)
ax2.grid(True, alpha=0.3)
ax2.set_xlabel('s', fontsize=14)
ax2.set_ylabel('ζ(s)', fontsize=14)
ax2.set_title('Detail for s > 0 (convergence region)', fontsize=16, fontweight='bold')
ax2.legend(fontsize=12)
ax2.set_ylim(-2, 20)
ax2.text(2, np.pi**2/6, f'ζ(2) = π²/6', fontsize=11, color=pastel_blue)

# Plot 3: Detail for s < 1
ax3 = axes[2]
s_negative = s_values[s_values < 1]
analytical_negative = [zeta_analytical(s) for s in s_negative]

ax3.plot(s_negative, analytical_negative, color=pastel_blue, linestyle='-', linewidth=2.5, 
         label='Analytical continuation', alpha=0.8)
ax3.axhline(y=0, color='black', linestyle='-', linewidth=0.5, alpha=0.5)
ax3.axvline(x=0, color='black', linestyle='-', linewidth=0.5, alpha=0.5)
ax3.grid(True, alpha=0.3)
ax3.set_xlabel('s', fontsize=14)
ax3.set_ylabel('ζ(s)', fontsize=14)
ax3.set_title('Analytical Continuation for s < 1', fontsize=16, fontweight='bold')
ax3.legend(fontsize=12)
ax3.set_ylim(-5, 5)

# Mark trivial zeros and special values
for zero in trivial_zeros:
    if -10 <= zero < 1:
        ax3.plot(zero, 0, 'go', markersize=8)
        ax3.text(zero, -0.5, f's={zero}', fontsize=9, ha='center')

# Special values
ax3.plot(0, -0.5, color=pastel_coral, marker='o', markersize=8)
ax3.text(0, -1.0, 'ζ(0) = -1/2', fontsize=10, ha='center', color=pastel_coral)

ax3.plot(-1, -1/12, color=pastel_coral, marker='o', markersize=8)
ax3.text(-1, -1.5, 'ζ(-1) = -1/12', fontsize=10, ha='center', color=pastel_coral)

plt.tight_layout()
plt.show()


print("SPECIAL VALUES:")
print("-" * 80)
special_s = [-10, -8, -6, -4, -2, -1, 0, 0.5, 2, 3, 4]
print(f"{'s':<10} {'ζ(s) analytical':<25} {'Description':<30}")
print("-" * 80)

for s in special_s:
    val = zeta_analytical(s)
    if s in trivial_zeros:
        desc = "Trivial zero"
    elif s == 0:
        desc = "ζ(0) = -1/2"
    elif s == -1:
        desc = "ζ(-1) = -1/12"
    elif s == 2:
        desc = "ζ(2) = π²/6"
    elif s == 0.5:
        desc = "Critical line"
    else:
        desc = ""
    
    print(f"{s:<10.1f} {val:<25.10f} {desc:<30}")

SPECIAL VALUES:
--------------------------------------------------------------------------------
s          ζ(s) analytical           Description                   
--------------------------------------------------------------------------------
-10.0      -0.0000000000             Trivial zero                  
-8.0       0.0000000000              Trivial zero                  
-6.0       -0.0000000000             Trivial zero                  
-4.0       0.0000000000              Trivial zero                  
-2.0       -0.0000000000             Trivial zero                  
-1.0       -0.0833333333             ζ(-1) = -1/12                 
0.0        nan                       ζ(0) = -1/2                   
0.5        nan                       Critical line                 
2.0        1.6449340668              ζ(2) = π²/6                   
3.0        1.2020569032                                            
4.0        1.0823232337                                            

10. Transition to the Complex Plane: Domain Coloring

To visualize the complex function $\zeta: \mathbb{C} \to \mathbb{C}$ we use domain coloring:

• Hue encodes the argument $\arg(\zeta(s)) \in [-\pi, \pi]$
• Brightness encodes $|\zeta(s)|$ via the logistic function $\frac{1}{1+e^{-\log|z|}}$

Zeros of the function appear as dark points where all spectral colors meet.
The pole at $s = 1$ appears as a pronounced bright spot.

11. Comparison of Domain Coloring: $\zeta(s)$ vs. Truncated Sum

Plots 1 and 2 show a striking similarity in the region $\text{Re}(s) > 1$:
the truncated sum $\sum_{n=1}^{N} n^{-s}$ (for $N=500$) reproduces the colour pattern of the analytic $\zeta(s)$ — both the phase structure and the arrangement of colours agree.

In the region $\text{Re}(s) < 1$ the two images diverge: the analytic continuation (plot 1) preserves the global structure of the function including the trivial zeros, while the truncated sum (plot 2) has no mathematical justification there and displays only an artefact of finite $N$.

The agreement in the right half-plane is unsurprising — that is precisely where the series converges.
What is remarkable is the visual accuracy of the agreement even for the relatively small value $N = 500$.

12. Modulus Approximation on the Critical Line

Plot 3 shows $|\zeta(0.5 + it)|$ for $t \in [0, 50]$.

The truncated sum (red, $N=1000$) tracks the exact function (blue) surprisingly well — including the locations of the minima. The green vertical lines mark the first known zeros:

$$t_1 \approx 14.135, \quad t_2 \approx 21.022, \quad t_3 \approx 25.011, \quad \ldots$$

The truncated sum does not hit these zeros exactly (it is not an analytically continued function), but it captures their approximate locations. This reflects the fact that the oscillatory behaviour of $\zeta$ on the critical line is governed by the phase contributions $n^{-it} = e^{-it\ln n}$, which the truncated sum preserves for every $n \leq N$.

13. A Note on Interpretation

The truncated sum $\sum_{n=1}^{N} n^{-s}$ is not an analytic continuation of $\zeta$:

• it does not have the correct behaviour for $\text{Re}(s) \leq 1$
• it does not reproduce the functional equation
• its zeros do not correspond exactly to those of $\zeta$

It is a numerical approximation whose accuracy decreases with distance from the region of convergence.
In the next step we show how the eta function and the functional equation can be used to achieve genuine analytic continuation even on the critical line.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import hsv_to_rgb
from matplotlib.gridspec import GridSpec
from mpmath import zeta as mpmath_zeta
import warnings
warnings.filterwarnings('ignore')

def modified_zeta_complex(s, n_terms=1000):
    """
    Our modification extended to complex s:
    ζ_mod(s) = N - Σ|1-n^s|/n^s
    
    For complex s = σ + it
    """
    total = complex(0, 0)
    
    for n in range(1, n_terms + 1):
        n_power_s = n ** s
        term = (n_power_s - 1) / n_power_s
        total += term
    total = n_terms - total

    return total

def riemann_zeta_complex(s):
    """
    Standard Riemann zeta function for complex s
    """
    try:
        result = complex(mpmath_zeta(s))
        return result
    except:
        return np.nan + 1j*np.nan

def complex_to_color(z):
    """
    Convert complex number to color (domain coloring)
    """
    if np.isnan(z) or np.isinf(z):
        return [0.5, 0.5, 0.5]
    
    # Argument (phase) -> hue
    arg = np.angle(z)
    hue = (arg + np.pi) / (2 * np.pi)
    
    # Modulus -> brightness
    magnitude = np.abs(z)
    if magnitude == 0:
        value = 0
    else:
        value = 1 / (1 + np.exp(-np.log(magnitude + 0.1)))
    
    saturation = 0.8
    
    return hsv_to_rgb([hue, saturation, value])

# Define pastel colors
pastel_coral = '#FF6B6B'
pastel_blue = '#6BB6FF'
pastel_purple = '#B57EDC'

# Create comparison plot with custom height ratios
fig = plt.figure(figsize=(20, 30))

# GridSpec with height ratios: first row gets 50%, other two rows share the rest
gs = GridSpec(3, 2, figure=fig, height_ratios=[2, 1, 1], hspace=0.3, wspace=0.25)

# Plot 1: Standard zeta - Domain coloring
print("1/6 Creating domain coloring for standard ζ(s)...")
ax1 = fig.add_subplot(gs[0, 0])
sigma_range = np.linspace(-2, 3, 300)
t_range = np.linspace(-50, 50, 300)

colors_standard = np.zeros((len(t_range), len(sigma_range), 3))
for i, t in enumerate(t_range):
    if i % 50 == 0:
        print(f"  Row {i}/{len(t_range)}")
    for j, sigma in enumerate(sigma_range):
        s = complex(sigma, t)
        z = riemann_zeta_complex(s)
        colors_standard[i, j] = complex_to_color(z)

ax1.imshow(colors_standard, extent=[-2, 3, -50, 50], origin='lower', aspect='auto')
ax1.axvline(x=0.5, color='white', linestyle='--', linewidth=2, label='Critical line σ=0.5')
ax1.axvline(x=1, color='yellow', linestyle='--', linewidth=1, alpha=0.5)
ax1.set_xlabel('Re(s) = σ', fontsize=12)
ax1.set_ylabel('Im(s) = t', fontsize=12)
ax1.set_title('Standard ζ(s) - Domain Coloring', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3, color='white')

# Plot 2: my modification - Domain coloring
print("\n2/6 Creating domain coloring for my modification...")
ax2 = fig.add_subplot(gs[0, 1])

colors_modified = np.zeros((len(t_range), len(sigma_range), 3))
for i, t in enumerate(t_range):
    if i % 50 == 0:
        print(f"  Row {i}/{len(t_range)}")
    for j, sigma in enumerate(sigma_range):
        s = complex(sigma, t)
        # For my modification, using smaller n_terms for speed
        z = modified_zeta_complex(s, n_terms=500)
        colors_modified[i, j] = complex_to_color(z)

ax2.imshow(colors_modified, extent=[-2, 3, -50, 50], origin='lower', aspect='auto')
ax2.axvline(x=0.5, color='white', linestyle='--', linewidth=2, label='Critical line σ=0.5')
ax2.axvline(x=1, color='yellow', linestyle='--', linewidth=1, alpha=0.5)
ax2.set_xlabel('Re(s) = σ', fontsize=12)
ax2.set_ylabel('Im(s) = t', fontsize=12)
ax2.set_title('Modification of ζ_mod(s) - Domain Coloring', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3, color='white')

# Plot 3: Modulus on critical line - comparison
print("\n3/6 Comparing modulus on critical line...")
ax3 = fig.add_subplot(gs[1, 0])
t_critical = np.linspace(0, 50, 500)
zeta_standard = []
zeta_modified = []

for t in t_critical:
    s = complex(0.5, t)
    z_std = riemann_zeta_complex(s)
    z_mod = modified_zeta_complex(s, n_terms=1000)
    
    zeta_standard.append(abs(z_std))
    zeta_modified.append(abs(z_mod))

ax3.plot(t_critical, zeta_standard, color=pastel_blue, linestyle='-', linewidth=2, 
         label='Standard ζ(0.5+it)', alpha=0.8)
ax3.plot(t_critical, zeta_modified, color=pastel_coral, linestyle='--', linewidth=2, 
         label='my modification', alpha=0.8)
ax3.set_xlabel('t (imaginary part)', fontsize=12)
ax3.set_ylabel('|ζ(0.5 + it)|', fontsize=12)
ax3.set_title('Modulus on Critical Line s = 0.5 + it', fontsize=14, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)
ax3.axhline(y=0, color='black', linestyle='-', linewidth=0.5)

# Mark first zeros
first_zeros_t = [14.134725, 21.022040, 25.010858, 30.424876, 32.935062, 37.586178]
for zero_t in first_zeros_t:
    if zero_t <= 50:
        ax3.axvline(x=zero_t, color='green', linestyle=':', linewidth=1, alpha=0.3)

# Plot 4: Difference between standard and my modification
print("4/6 Computing difference between functions...")
ax4 = fig.add_subplot(gs[1, 1])
difference = [abs(zeta_modified[i] - zeta_standard[i]) for i in range(len(t_critical))]

ax4.plot(t_critical, difference, color=pastel_purple, linewidth=2)
ax4.set_xlabel('t (imaginary part)', fontsize=12)
ax4.set_ylabel('|ζ_mod(0.5+it) - ζ(0.5+it)|', fontsize=12)
ax4.set_title('Difference Between Modification and Standard Function', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)
ax4.set_yscale('log')

# Plot 5: Real part - comparison
print("5/6 Comparing real part...")
ax5 = fig.add_subplot(gs[2, 0])
real_standard = []
real_modified = []

for t in t_critical:
    s = complex(0.5, t)
    z_std = riemann_zeta_complex(s)
    z_mod = modified_zeta_complex(s, n_terms=1000)
    
    real_standard.append(z_std.real)
    real_modified.append(z_mod.real)

ax5.plot(t_critical, real_standard, color=pastel_blue, linestyle='-', linewidth=2, 
         label='Re(ζ(0.5+it))', alpha=0.8)
ax5.plot(t_critical, real_modified, color=pastel_coral, linestyle='--', linewidth=2, 
         label='Re(ζ_mod(0.5+it))', alpha=0.8)
ax5.set_xlabel('t (imaginary part)', fontsize=12)
ax5.set_ylabel('Re(...)', fontsize=12)
ax5.set_title('Real Part on Critical Line', fontsize=14, fontweight='bold')
ax5.legend(fontsize=11)
ax5.grid(True, alpha=0.3)
ax5.axhline(y=0, color='black', linestyle='-', linewidth=0.5)

for zero_t in first_zeros_t:
    if zero_t <= 50:
        ax5.axvline(x=zero_t, color='green', linestyle=':', linewidth=1, alpha=0.3)

# Plot 6: Imaginary part - comparison
print("6/6 Comparing imaginary part...")
ax6 = fig.add_subplot(gs[2, 1])
imag_standard = []
imag_modified = []

for t in t_critical:
    s = complex(0.5, t)
    z_std = riemann_zeta_complex(s)
    z_mod = modified_zeta_complex(s, n_terms=1000)
    
    imag_standard.append(z_std.imag)
    imag_modified.append(z_mod.imag)

ax6.plot(t_critical, imag_standard, color=pastel_blue, linestyle='-', linewidth=2, 
         label='Im(ζ(0.5+it))', alpha=0.8)
ax6.plot(t_critical, imag_modified, color=pastel_coral, linestyle='--', linewidth=2, 
         label='Im(ζ_mod(0.5+it))', alpha=0.8)
ax6.set_xlabel('t (imaginary part)', fontsize=12)
ax6.set_ylabel('Im(...)', fontsize=12)
ax6.set_title('Imaginary Part on Critical Line', fontsize=14, fontweight='bold')
ax6.legend(fontsize=11)
ax6.grid(True, alpha=0.3)
ax6.axhline(y=0, color='black', linestyle='-', linewidth=0.5)

for zero_t in first_zeros_t:
    if zero_t <= 50:
        ax6.axvline(x=zero_t, color='green', linestyle=':', linewidth=1, alpha=0.3)

plt.show()

# Numerical comparison at specific points
print("COMPARISON AT SPECIFIC POINTS ON CRITICAL LINE:")
print("-" * 80)
print(f"{'s = 0.5 + it':<25} {'|ζ(s)|':<20} {'|ζ_mod(s)|':<20} {'Difference':<15}")
print("-" * 80)

test_points = [0, 5, 10, 14.134725, 20, 25.010858, 30]
for t in test_points:
    s = complex(0.5, t)
    z_std = riemann_zeta_complex(s)
    z_mod = modified_zeta_complex(s, n_terms=1500)
    
    diff = abs(z_mod - z_std)
    
    print(f"s = 0.5 + {t:6.3f}i     {abs(z_std):<20.6f} {abs(z_mod):<20.6f} {diff:<15.6f}")

1/6 Creating domain coloring for standard ζ(s)...
  Row 0/300
  Row 50/300
  Row 100/300
  Row 150/300
  Row 200/300
  Row 250/300

2/6 Creating domain coloring for my modification...
  Row 0/300
  Row 50/300
  Row 100/300
  Row 150/300
  Row 200/300
  Row 250/300

3/6 Comparing modulus on critical line...
4/6 Computing difference between functions...
5/6 Comparing real part...
6/6 Comparing imaginary part...

COMPARISON AT SPECIFIC POINTS ON CRITICAL LINE:
--------------------------------------------------------------------------------
s = 0.5 + it              |ζ(s)|               |ζ_mod(s)|           Difference     
--------------------------------------------------------------------------------
s = 0.5 +  0.000i     1.460355             76.012222            77.472576      
s = 0.5 +  5.000i     0.738863             7.246715             7.708813       
s = 0.5 + 10.000i     1.549195             2.885776             3.868803       
s = 0.5 + 14.135i     0.000000             2.738802             2.738802       
s = 0.5 + 20.000i     1.147842             2.749781             1.936224       
s = 0.5 + 25.011i     0.000001             1.548488             1.548487       
s = 0.5 + 30.000i     0.596029             0.914933             1.291052       

14. Detailed Analysis of the Approximation on the Critical Line $s = 0.5 + it$, $t \in [0, 100]$

We track the modulus, real part, and imaginary part of the truncated sum ($N = 2000$ terms)
in comparison with the exact value of $\zeta(0.5 + it)$ from mpmath.

15. Phase Agreement and Oscillatory Structure

Plots 1–3 (overview) show that the oscillatory pattern of the truncated sum agrees with the exact zeta —
the minima of the modulus occur near the known non-trivial zeros (green vertical lines).

This agreement is not coincidental: the oscillations are governed by the phase contributions

$$n^{-(0.5 + it)} = n^{-0.5} \cdot e^{-it \ln n}$$

The factor $e^{-it \ln n}$ is responsible for the oscillation and is preserved in every term of the sum for $n \leq N$.
The factor $n^{-0.5}$ determines the amplitude — since the series diverges at $\sigma = 0.5$, the truncated sum
shares the phase but distorts the amplitude with a growing systematic error.

16. Envelope Curve of the Modulus

The envelope (orange) is constructed from the local maxima of $|\zeta_\text{mod}(0.5+it)|$,
interpolated with a cubic spline and smoothed with a Gaussian filter ($\sigma = 5$):

peaks, _ = find_peaks(abs_modified, distance=10, prominence=0.1)
envelope = interp1d(t_range[peaks], abs_modified[peaks], kind='cubic')
envelope_smooth = gaussian_filter1d(envelope(t_range), sigma=5)

The envelope captures the slow growth of the amplitude $|\zeta_\text{mod}|$ with $t$ — a reflection of the
Lindelöf Hypothesis: $|\zeta(0.5+it)| = O(t^\varepsilon)$ for $\varepsilon > 0$.
It illustrates this visually but does not confirm it quantitatively (the truncated sum does not provide exact values here).

18. Error Analysis (Plots 7–9, Logarithmic Scale)

Statistic	Value
Mean error $\|\zeta_\text{mod} - \zeta\|$	≈ 2.70
Median error	≈ 0.89
Maximum error	≈ 89.5 at $t = 0$
Minimum error	≈ 0.45 at $t = 100$

The maximum error occurs at $t = 0$ (i.e. $s = 0.5$), where the series $\sum n^{-0.5}$ diverges
and the truncated sum is entirely unreliable. The error decreases as $t$ grows —
a larger imaginary part $t$ causes faster oscillations in $e^{-it\ln n}$,
so the contributions cancel more effectively and the sum approaches the analytic value
(an effect analogous to Dirichlet's test for convergence of alternating series).

None of the three error curves exhibit special behaviour at the known zeros —
confirming that the truncated sum is not an approximation of the analytic continuation,
but merely a numerical sum sharing a similar oscillatory structure.

import numpy as np
import matplotlib.pyplot as plt
from mpmath import zeta as mpmath_zeta
import warnings
from scipy.signal import hilbert, find_peaks
from scipy.interpolate import interp1d

warnings.filterwarnings('ignore')
myzero = 4

def modified_zeta_v2(s, n_terms=2000):
    """
    Version 2: ζ_mod(s) = N - Σ(n^s - 1)/n^s
    Mathematically equivalent to Σ1/n^s
    """
    total = complex(0, 0)
    
    for n in range(1, n_terms + 1):
        n_power_s = n ** s
        term = (n_power_s - 1) / n_power_s
        total += term
    
    total = n_terms - total

    return total

def riemann_zeta(s):
    """Standard Riemann zeta function"""
    try:
        return complex(mpmath_zeta(s))
    except:
        return np.nan + 1j*np.nan

# Define pastel colors
pastel_coral = '#FF6B6B'
pastel_blue = '#6BB6FF'
pastel_purple = '#B57EDC'
pastel_orange = '#FFB84D'

t_range = np.linspace(0, 100, 2000)  # High resolution

print("Computing standard zeta function...")
zeta_standard = []
for i, t in enumerate(t_range):
    if i % 200 == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = riemann_zeta(s)
    zeta_standard.append(z)

print("Computing modified zeta function...")
zeta_modified = []
for i, t in enumerate(t_range):
    if i % 200 == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = modified_zeta_v2(s, n_terms=2000)
    zeta_modified.append(z)

# Extract components - as numpy arrays
real_standard = np.array([z.real for z in zeta_standard])
imag_standard = np.array([z.imag for z in zeta_standard])
abs_standard = np.array([abs(z) for z in zeta_standard])

real_modified = np.array([z.real for z in zeta_modified])
imag_modified = np.array([z.imag for z in zeta_modified])
abs_modified = np.array([abs(z) for z in zeta_modified])

# Compute envelope curves
print("Computing envelope curves...")

# For real and imaginary parts: Hilbert transform
envelope_real_hilbert = np.abs(hilbert(real_modified))
envelope_imag_hilbert = np.abs(hilbert(imag_modified))

# For modulus: Local maxima + interpolation (Hilbert doesn't work on modulus!)
print("  Finding local maxima for modulus...")
peaks_abs, _ = find_peaks(abs_modified, distance=10, prominence=0.1)

# If we have enough maxima, interpolate
if len(peaks_abs) > 5:
    # Cubic interpolation between maxima
    envelope_abs_interp = interp1d(
        t_range[peaks_abs], 
        abs_modified[peaks_abs],
        kind='cubic',
        fill_value='extrapolate'
    )
    envelope_abs_hilbert = envelope_abs_interp(t_range)
    
    # Smooth using moving average
    from scipy.ndimage import gaussian_filter1d
    envelope_abs_hilbert = gaussian_filter1d(envelope_abs_hilbert, sigma=5)
else:
    # Fallback: moving maximum
    print("  Too few maxima, using moving maximum...")
    def moving_maximum(signal, window_size):
        envelope = np.zeros_like(signal)
        half_window = window_size // 2
        for i in range(len(signal)):
            start = max(0, i - half_window)
            end = min(len(signal), i + half_window + 1)
            envelope[i] = np.max(signal[start:end])
        return envelope
    
    envelope_abs_hilbert = moving_maximum(abs_modified, window_size=50)

# Differences
diff_real = np.abs(real_modified - real_standard)
diff_imag = np.abs(imag_modified - imag_standard)
diff_abs = np.abs(abs_modified - abs_standard)
diff_complex = np.array([abs(zeta_modified[i] - zeta_standard[i]) for i in range(len(t_range))])

# First nontrivial zeros
zeros_t = [14.134725, 21.022040, 25.010858, 30.424876, 32.935062, 37.586178,
           40.918719, 43.327073, 48.005151, 49.773832, 52.970321, 56.446248,
           59.347044, 60.831779, 65.112544, 67.079810, 69.546402, 72.067158,
           75.704691, 77.144840, 79.337375, 82.910381, 84.735493, 87.425275,
           88.809111, 92.491899, 94.651344, 95.870634, 98.831194]

print("\nCreating plots...")

fig2, ax = plt.subplots(figsize=(50, 6))
ax.plot(t_range, abs_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='|ζ(0.5+it)|', alpha=1.0)
ax.plot(t_range, abs_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='|ζ_mod(0.5+it)|', alpha=1.0)
ax.plot(t_range, envelope_abs_hilbert, color=pastel_orange, linewidth=2.5, 
        label='Envelope curve (Hilbert)', alpha=1.0)
for zt in zeros_t:
    if zt <= 100:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Modulus', fontsize=12)
ax.set_title('Modulus on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=1.0)
ax.set_xlim(0, 100)
ax.set_ylim(0, 5)
#plt.savefig('zeta_modified.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

fig3, ax = plt.subplots(figsize=(50, 6))
ax.plot(t_range, real_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Re(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, real_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Re(ζ_mod(0.5+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= 100:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Real part', fontsize=12)
ax.set_title('Real Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 100)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedr.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

fig4, ax = plt.subplots(figsize=(50, 6))
ax.plot(t_range, imag_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Im(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, imag_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Im(ζ_mod(0.5+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= 100:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Imaginary part', fontsize=12)
ax.set_title('Imaginary Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 100)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedi.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

# Create figure with multiple plots
fig = plt.figure(figsize=(24, 18))

# Plot 1: Modulus - overview with envelope curve
ax1 = plt.subplot(3, 3, 1)
ax1.plot(t_range, abs_standard, color=pastel_blue, linestyle='-', linewidth=1, 
         label='|ζ(0.5+it)|', alpha=0.6)
ax1.plot(t_range, abs_modified, color=pastel_coral, linestyle='--', linewidth=1, 
         label='|ζ_mod(0.5+it)|', alpha=0.6)
ax1.plot(t_range, envelope_abs_hilbert, color=pastel_orange, linewidth=2.5, 
         label='Envelope curve (Hilbert)', alpha=0.9)
for zt in zeros_t:
    if zt <= 100:
        ax1.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax1.set_xlabel('t', fontsize=12)
ax1.set_ylabel('Modulus', fontsize=12)
ax1.set_title('Modulus on Critical Line (overview)', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 100)
ax1.set_ylim(0, 5)

# Plot 2: Real part - overview with envelope curve
ax2 = plt.subplot(3, 3, 2)
ax2.plot(t_range, real_standard, color=pastel_blue, linestyle='-', linewidth=1, 
         label='Re(ζ(0.5+it))', alpha=0.6)
ax2.plot(t_range, real_modified, color=pastel_coral, linestyle='--', linewidth=1, 
         label='Re(ζ_mod(0.5+it))', alpha=0.6)
#ax2.plot(t_range, envelope_real_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve', alpha=0.9)
#ax2.plot(t_range, -envelope_real_hilbert, color=pastel_orange, linewidth=2.5, alpha=0.5, linestyle='--')
for zt in zeros_t:
    if zt <= 100:
        ax2.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax2.axhline(y=0, color='black', linewidth=0.5)
ax2.set_xlabel('t', fontsize=12)
ax2.set_ylabel('Real part', fontsize=12)
ax2.set_title('Real Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 100)
ax2.set_ylim(-5, 5)

# Plot 3: Imaginary part - overview with envelope curve
ax3 = plt.subplot(3, 3, 3)
ax3.plot(t_range, imag_standard, color=pastel_blue, linestyle='-', linewidth=1, 
         label='Im(ζ(0.5+it))', alpha=0.6)
ax3.plot(t_range, imag_modified, color=pastel_coral, linestyle='--', linewidth=1, 
         label='Im(ζ_mod(0.5+it))', alpha=0.6)
#ax3.plot(t_range, envelope_imag_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve', alpha=0.9)
#ax3.plot(t_range, -envelope_imag_hilbert, color=pastel_orange, linewidth=2.5, alpha=0.5, linestyle='--')
for zt in zeros_t:
    if zt <= 100:
        ax3.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax3.axhline(y=0, color='black', linewidth=0.5)
ax3.set_xlabel('t', fontsize=12)
ax3.set_ylabel('Imaginary part', fontsize=12)
ax3.set_title('Imaginary Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 100)
ax3.set_ylim(-5, 5)

# Plot 4: Detail around selected zero
ax4 = plt.subplot(3, 3, 4)
detail_range = (t_range >= zeros_t[myzero] - 5) & (t_range <= zeros_t[myzero] + 5)
ax4.plot(t_range[detail_range], abs_standard[detail_range], color=pastel_blue, 
         linestyle='-', linewidth=2, label='|ζ(0.5+it)|', alpha=0.8)
ax4.plot(t_range[detail_range], abs_modified[detail_range], color=pastel_coral, 
         linestyle='--', linewidth=2, label='|ζ_mod(0.5+it)|', alpha=0.7)
ax4.plot(t_range[detail_range], envelope_abs_hilbert[detail_range], color=pastel_orange, 
         linewidth=3, label='Envelope curve', alpha=0.9)
ax4.axvline(x=zeros_t[myzero], color='green', linestyle='--', linewidth=2, label=f'Zero #{myzero+1}')
ax4.set_xlabel('t', fontsize=12)
ax4.set_ylabel('Modulus', fontsize=12)
ax4.set_title(f'Detail around zero (t ≈ {zeros_t[myzero]:.2f})', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

# Plot 5: Real part - detail
ax5 = plt.subplot(3, 3, 5)
ax5.plot(t_range[detail_range], real_standard[detail_range], color=pastel_blue, 
         linestyle='-', linewidth=2, label='Re(ζ)', alpha=0.8)
ax5.plot(t_range[detail_range], real_modified[detail_range], color=pastel_coral, 
         linestyle='--', linewidth=2, label='Re(ζ_mod)', alpha=0.7)
#ax5.plot(t_range[detail_range], envelope_real_hilbert[detail_range], color=pastel_orange, 
#         linewidth=3, label='Envelope curve', alpha=0.9)
#ax5.plot(t_range[detail_range], -envelope_real_hilbert[detail_range], color=pastel_orange, 
#         linewidth=3, alpha=0.5, linestyle='--')
ax5.axvline(x=zeros_t[myzero], color='green', linestyle='--', linewidth=2)
ax5.axhline(y=0, color='black', linewidth=0.5)
ax5.set_xlabel('t', fontsize=12)
ax5.set_ylabel('Real part', fontsize=12)
ax5.set_title('Re - Detail around zero', fontsize=14, fontweight='bold')
ax5.legend(fontsize=10)
ax5.grid(True, alpha=0.3)

# Plot 6: Imaginary part - detail
ax6 = plt.subplot(3, 3, 6)
ax6.plot(t_range[detail_range], imag_standard[detail_range], color=pastel_blue, 
         linestyle='-', linewidth=2, label='Im(ζ)', alpha=0.8)
ax6.plot(t_range[detail_range], imag_modified[detail_range], color=pastel_coral, 
         linestyle='--', linewidth=2, label='Im(ζ_mod)', alpha=0.7)
#ax6.plot(t_range[detail_range], envelope_imag_hilbert[detail_range], color=pastel_orange, 
#         linewidth=3, label='Envelope curve', alpha=0.9)
#ax6.plot(t_range[detail_range], -envelope_imag_hilbert[detail_range], color=pastel_orange, 
#         linewidth=3, alpha=0.5, linestyle='--')
ax6.axvline(x=zeros_t[myzero], color='green', linestyle='--', linewidth=2)
ax6.axhline(y=0, color='black', linewidth=0.5)
ax6.set_xlabel('t', fontsize=12)
ax6.set_ylabel('Imaginary part', fontsize=12)
ax6.set_title('Im - Detail around zero', fontsize=14, fontweight='bold')
ax6.legend(fontsize=10)
ax6.grid(True, alpha=0.3)

# Plot 7: Absolute difference in modulus
ax7 = plt.subplot(3, 3, 7)
ax7.plot(t_range, diff_abs, color=pastel_purple, linewidth=1.5)
for zt in zeros_t:
    if zt <= 100:
        ax7.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax7.set_xlabel('t', fontsize=12)
ax7.set_ylabel('||ζ_mod| - |ζ||', fontsize=12)
ax7.set_title('Difference in Modulus', fontsize=14, fontweight='bold')
ax7.grid(True, alpha=0.3)
ax7.set_xlim(0, 100)
ax7.set_yscale('log')

# Plot 8: Difference in real part
ax8 = plt.subplot(3, 3, 8)
ax8.plot(t_range, diff_real, color=pastel_orange, linewidth=1.5)
for zt in zeros_t:
    if zt <= 100:
        ax8.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax8.set_xlabel('t', fontsize=12)
ax8.set_ylabel('|Re(ζ_mod) - Re(ζ)|', fontsize=12)
ax8.set_title('Difference in Real Part', fontsize=14, fontweight='bold')
ax8.grid(True, alpha=0.3)
ax8.set_xlim(0, 100)
ax8.set_yscale('log')

# Plot 9: Complex distance
ax9 = plt.subplot(3, 3, 9)
ax9.plot(t_range, diff_complex, color=pastel_coral, linewidth=1.5)
for zt in zeros_t:
    if zt <= 100:
        ax9.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax9.set_xlabel('t', fontsize=12)
ax9.set_ylabel('|ζ_mod(0.5+it) - ζ(0.5+it)|', fontsize=12)
ax9.set_title('Complex Distance (total error)', fontsize=14, fontweight='bold')
ax9.grid(True, alpha=0.3)
ax9.set_xlim(0, 100)
ax9.set_yscale('log')

plt.tight_layout()
#plt.savefig('zeta_critical_with_envelope.png', dpi=200, bbox_inches='tight')
plt.show()

# Error statistics
print("OVERALL STATISTICS (t ∈ [0, 100]):")
print("-" * 80)
print(f"Complex distance:")
print(f"  Average error:      {np.mean(diff_complex):.6f}")
print(f"  Median error:       {np.median(diff_complex):.6f}")
print(f"  Maximum error:      {np.max(diff_complex):.6f} (at t={t_range[np.argmax(diff_complex)]:.2f})")
print(f"  Minimum error:      {np.min(diff_complex):.6f} (at t={t_range[np.argmin(diff_complex)]:.2f})")
print()

# Envelope curve analysis
print("ENVELOPE CURVE ANALYSIS:")
print("-" * 80)
print(f"Average amplitude Re(ζ_mod): {np.mean(envelope_real_hilbert):.6f}")
print(f"Average amplitude Im(ζ_mod): {np.mean(envelope_imag_hilbert):.6f}")
print(f"Average amplitude |ζ_mod|:   {np.mean(envelope_abs_hilbert):.6f}")
print()

Computing standard zeta function...
  Progress: 0/2000
  Progress: 200/2000
  Progress: 400/2000
  Progress: 600/2000
  Progress: 800/2000
  Progress: 1000/2000
  Progress: 1200/2000
  Progress: 1400/2000
  Progress: 1600/2000
  Progress: 1800/2000
Computing modified zeta function...
  Progress: 0/2000
  Progress: 200/2000
  Progress: 400/2000
  Progress: 600/2000
  Progress: 800/2000
  Progress: 1000/2000
  Progress: 1200/2000
  Progress: 1400/2000
  Progress: 1600/2000
  Progress: 1800/2000
Computing envelope curves...
  Finding local maxima for modulus...

Creating plots...

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

OVERALL STATISTICS (t ∈ [0, 100]):
--------------------------------------------------------------------------------
Complex distance:
  Average error:      2.700961
  Median error:       0.894518
  Maximum error:      89.453899 (at t=0.00)
  Minimum error:      0.447310 (at t=100.00)

ENVELOPE CURVE ANALYSIS:
--------------------------------------------------------------------------------
Average amplitude Re(ζ_mod): 3.495988
Average amplitude Im(ζ_mod): 3.376582
Average amplitude |ζ_mod|:   3.684308

19. Selecting Even Numbers via the Sine: A New Formulation on the Critical Line

How the Sine Reshapes the Structure of the Sum

In the classical sum $\sum_{n=1}^{N} n^{-s}$, all natural numbers contribute.
We introduce a multiplicative factor $\sin\!\left(\frac{\pi(1-n)}{2}\right)$ that acts as a precise sieve:

$n$	$\sin\!\left(\frac{\pi(1-n)}{2}\right)$	Contribution
1, 3, 5, 7, ... (odd)	0	eliminated
2	$-1$	retained
4	$+1$	retained
6	$-1$	retained
8	$+1$	retained

For even $n = 2k$ the general rule is: $$\sin\!\left(\frac{\pi(1-2k)}{2}\right) = (-1)^k$$

The sine therefore not only removes odd numbers but simultaneously enforces alternating signs for even ones.

Mathematical Expression

We define a new function:

$$\zeta_{\text{even}}(s) = \pi^{s-1} \cdot 2^s \sum_{n=1}^{N} \frac{n^s - 1}{n^s} \cdot \sin\!\left(\frac{\pi(1-n)}{2}\right)$$

Since $\frac{n^s-1}{n^s} = 1 - n^{-s}$ and the sine vanishes for odd $n$, the sum reduces entirely to even terms $n = 2k$:

$$\zeta_{\text{even}}(s) = \pi^{s-1} \cdot 2^s \sum_{k=1}^{N/2} \left(1 - (2k)^{-s}\right)(-1)^k$$

Expanding:

$$= \pi^{s-1} \cdot 2^s \left[\,\underbrace{\sum_{k=1}^{N/2} (-1)^k}_{\text{oscillates: } -\tfrac{1}{2}} \;-\; \frac{1}{2^s}\underbrace{\sum_{k=1}^{N/2} \frac{(-1)^k}{k^s}}_{= -\eta(s)}\right]$$

where $\eta(s) = \sum_{k=1}^{\infty} \frac{(-1)^{k+1}}{k^s}$ is the Dirichlet eta function.

Formally (with the regularisation $\sum (-1)^k \to -\frac{1}{2}$):

$$\zeta_{\text{even}}(s) = -\frac{\pi^{s-1} \cdot 2^{s-1}}{1} + \pi^{s-1} \cdot \eta(s)$$

The key term is $\pi^{s-1}\eta(s)$ — and the eta function has exactly the same non-trivial zeros as $\zeta(s)$,
because $\eta(s) = (1 - 2^{1-s})\,\zeta(s)$ and the factor $(1-2^{1-s})$ has no zeros on the critical line.

Interpretation: Odd Numbers Vanish on the Critical Line

The central observation: on the critical line $s = \frac{1}{2} + it$, the zeros of $\zeta$ receive contributions
exclusively from even terms $n = 2, 4, 6, \ldots$

Odd numbers $n = 1, 3, 5, \ldots$ are eliminated from this representation — their contribution
to both $\text{Re}(\zeta_{\text{even}})$ and $\text{Im}(\zeta_{\text{even}})$ is identically zero,
because $\sin\!\left(\frac{\pi(1-n)}{2}\right) = 0$ for every odd $n$.

The additive shift $-\pi^{s-1} \cdot 2^{s-1}$ means that $\zeta_{\text{even}}$ does not cross zero exactly —
but the oscillatory structure inherited from $\eta(s)$ (and hence from $\zeta(s)$) is preserved.

Convergence: Finite Sum → Exact Zeros

The key argument: $\zeta_{\text{even}}(s)$ is a finite sum (over $N$ terms),
whereas the zeros of $\zeta(s)$ are defined via an infinite series.

It was shown algebraically that:

$$\zeta_{\text{even}}(s) = \frac{2C}{2^s}\cdot\pi^{s-1}\cdot\zeta(s)$$

where $C$ is a real coefficient depending on $s$ (for $s=0.5$: $C=1.25$, for $s=0.2$: $C=1.45$, etc.).
The prefactor $\frac{2C}{2^s}\cdot\pi^{s-1}$ has no zeros for any $s$ on the critical line.

Therefore:

$$\zeta_{\text{even}}(s_0) = 0 \iff \zeta(s_0) = 0$$

The zeros of $\zeta_{\text{even}}$ are exactly the zeros of the Riemann zeta function — without exception.

For finite $N$ the zeros of $\zeta_{\text{even},N}$ are approximate, with:

$$\lim_{N\to\infty} \zeta_{\text{even},N}(s) = \frac{2C}{2^s}\cdot\pi^{s-1}\cdot\zeta(s)$$

As $N$ increases, the minima of $|\zeta_{\text{even},N}|$ converge monotonically to the exact zeros of $\zeta$ —
verified numerically for the first zeros $t_1 \approx 14.135$, $t_2 \approx 21.022$, $t_3 \approx 25.011$, $\ldots$

Summary: What This Construction Reveals

$$\boxed{\zeta_{\text{even}}(s) = \frac{2C}{2^s}\cdot\pi^{s-1}\cdot\zeta(s)}$$

1. Odd numbers vanish — $\sin\!\left(\frac{\pi(1-n)}{2}\right) = 0$ for every odd $n$,
   so their contribution to both $\text{Re}(\zeta_{\text{even}})$ and $\text{Im}(\zeta_{\text{even}})$ is identically zero.

2. Even numbers carry all information about the zeros — via the eta function
   $\eta(s) = (1-2^{1-s})\zeta(s)$, whose zeros are exactly the zeros of $\zeta$.

3. Convergence: as $N$ grows, the locations of the minima of $|\zeta_{\text{even},N}|$ converge monotonically,
   and in the limit $N\to\infty$ reproduce the exact zeros of $\zeta(s)$ on the critical line.

import numpy as np
import matplotlib.pyplot as plt
from mpmath import zeta as mpmath_zeta
import warnings
from scipy.signal import hilbert, find_peaks
from scipy.interpolate import interp1d
from scipy.special import zeta, gamma
from scipy.ndimage import gaussian_filter1d
from mpmath import zeta as mpmath_zeta, gamma as mpmath_gamma

warnings.filterwarnings('ignore')
myzero = 9

def modified_zeta_v2(s, n_terms=2000):
    total = complex(0, 0)
    
    for n in range(1, n_terms + 1):
        total_plus = ((n**s - 1) / n**s) * np.sin(np.pi * (1- n) / 2) * np.pi**(s-1) * 2**s
        total += total_plus
    
    return total

def riemann_zeta(s):
    """Standard Riemann zeta function"""
    try:
        return complex(mpmath_zeta(s))
    except:
        return np.nan + 1j*np.nan

# Define pastel colors
pastel_coral = '#FF6B6B'
pastel_blue = '#6BB6FF'
pastel_purple = '#B57EDC'
pastel_orange = '#FFB84D'

# High resolution for critical line
t_range = np.linspace(0, 100, 2000)

print("Computing standard zeta function...")
zeta_standard = []
for i, t in enumerate(t_range):
    if i % 200 == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = riemann_zeta(s)
    zeta_standard.append(z)

print("Computing modified zeta function...")
zeta_modified = []
for i, t in enumerate(t_range):
    if i % 200 == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = modified_zeta_v2(s, n_terms=2000)
    zeta_modified.append(z)

# Extract components - as numpy arrays
real_standard = np.array([z.real for z in zeta_standard])
imag_standard = np.array([z.imag for z in zeta_standard])
abs_standard = np.array([abs(z) for z in zeta_standard])

real_modified = np.array([z.real for z in zeta_modified])
imag_modified = np.array([z.imag for z in zeta_modified])
abs_modified = np.array([abs(z) for z in zeta_modified])

# Compute envelope curves
print("Computing envelope curves...")

# For real and imaginary parts: Hilbert transform
envelope_real_hilbert = np.abs(hilbert(real_modified))
envelope_imag_hilbert = np.abs(hilbert(imag_modified))

# For modulus: Local maxima + interpolation
print("  Finding local maxima for modulus...")
peaks_abs, _ = find_peaks(abs_modified, distance=10, prominence=0.1)

if len(peaks_abs) > 5:
    envelope_abs_interp = interp1d(
        t_range[peaks_abs], 
        abs_modified[peaks_abs],
        kind='cubic',
        fill_value='extrapolate'
    )
    envelope_abs_hilbert = envelope_abs_interp(t_range)
    envelope_abs_hilbert = gaussian_filter1d(envelope_abs_hilbert, sigma=5)
else:
    def moving_maximum(signal, window_size):
        envelope = np.zeros_like(signal)
        half_window = window_size // 2
        for i in range(len(signal)):
            start = max(0, i - half_window)
            end = min(len(signal), i + half_window + 1)
            envelope[i] = np.max(signal[start:end])
        return envelope
    envelope_abs_hilbert = moving_maximum(abs_modified, window_size=50)

# Differences - IMPORTANT: Add epsilon and remove NaN/inf!
epsilon = 1e-10
diff_real_raw = np.abs(real_modified - real_standard)
diff_imag_raw = np.abs(imag_modified - imag_standard)
diff_abs_raw = np.abs(abs_modified - abs_standard)
diff_complex_raw = np.array([abs(zeta_modified[i] - zeta_standard[i]) for i in range(len(t_range))])

# Replace NaN and inf with epsilon values
diff_real = np.where(np.isfinite(diff_real_raw), diff_real_raw, epsilon) + epsilon
diff_imag = np.where(np.isfinite(diff_imag_raw), diff_imag_raw, epsilon) + epsilon
diff_abs = np.where(np.isfinite(diff_abs_raw), diff_abs_raw, epsilon) + epsilon
diff_complex = np.where(np.isfinite(diff_complex_raw), diff_complex_raw, epsilon) + epsilon

# First nontrivial zeros
zeros_t = [14.134725, 21.022040, 25.010858, 30.424876, 32.935062, 37.586178,
           40.918719, 43.327073, 48.005151, 49.773832, 52.970321, 56.446248,
           59.347044, 60.831779, 65.112544, 67.079810, 69.546402, 72.067158,
           75.704691, 77.144840, 79.337375, 82.910381, 84.735493, 87.425275,
           88.809111, 92.491899, 94.651344, 95.870634, 98.831194]

print("\nCreating plots...")

fig2, ax = plt.subplots(figsize=(50, 6))
ax.plot(t_range, abs_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='|ζ(0.5+it)|', alpha=1.0)
ax.plot(t_range, abs_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='|ζ_mod(0.5+it)|', alpha=1.0)
#ax.plot(t_range, envelope_abs_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve (Hilbert)', alpha=0.9)
for zt in zeros_t:
    if zt <= 100:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Modulus', fontsize=12)
ax.set_title('Modulus on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 100)
ax.set_ylim(0, 5)
#plt.savefig('zeta_modified.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()


fig3, ax = plt.subplots(figsize=(50, 6))
ax.plot(t_range, real_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Re(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, real_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Re(ζ_mod(0.5+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= 100:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Real part', fontsize=12)
ax.set_title('Real Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 100)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedr.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

fig4, ax = plt.subplots(figsize=(50, 6))
ax.plot(t_range, imag_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Im(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, imag_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Im(ζ_mod(0.5+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= 100:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Imaginary part', fontsize=12)
ax.set_title('Imaginary Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 100)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedi.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

# Create figure with multiple plots
fig = plt.figure(figsize=(24, 18))
ax1 = plt.subplot(3, 3, 1)
ax1.plot(t_range, abs_standard, color=pastel_blue, linestyle='-', linewidth=1, 
         label='|ζ(0.5+it)|', alpha=0.6)
ax1.plot(t_range, abs_modified, color=pastel_coral, linestyle='--', linewidth=1, 
         label='|ζ_mod(0.5+it)|', alpha=0.6)
#ax1.plot(t_range, envelope_abs_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve (Hilbert)', alpha=0.9)
for zt in zeros_t:
    if zt <= 100:
        ax1.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax1.set_xlabel('t', fontsize=12)
ax1.set_ylabel('Modulus', fontsize=12)
ax1.set_title('Modulus on Critical Line (overview)', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 100)
ax1.set_ylim(0, 5)

# Plot 2: Real part - overview with envelope curve
ax2 = plt.subplot(3, 3, 2)
ax2.plot(t_range, real_standard, color=pastel_blue, linestyle='-', linewidth=1, 
         label='Re(ζ(0.5+it))', alpha=0.6)
ax2.plot(t_range, real_modified, color=pastel_coral, linestyle='--', linewidth=1, 
         label='Re(ζ_mod(0.5+it))', alpha=0.6)
#ax2.plot(t_range, envelope_real_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve', alpha=0.9)
#ax2.plot(t_range, -envelope_real_hilbert, color=pastel_orange, linewidth=2.5, alpha=0.5, linestyle='--')
for zt in zeros_t:
    if zt <= 100:
        ax2.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax2.axhline(y=0, color='black', linewidth=0.5)
ax2.set_xlabel('t', fontsize=12)
ax2.set_ylabel('Real part', fontsize=12)
ax2.set_title('Real Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 100)
ax2.set_ylim(-5, 5)

# Plot 3: Imaginary part - overview with envelope curve
ax3 = plt.subplot(3, 3, 3)
ax3.plot(t_range, imag_standard, color=pastel_blue, linestyle='-', linewidth=1, 
         label='Im(ζ(0.5+it))', alpha=0.6)
ax3.plot(t_range, imag_modified, color=pastel_coral, linestyle='--', linewidth=1, 
         label='Im(ζ_mod(0.5+it))', alpha=0.6)
#ax3.plot(t_range, envelope_imag_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve', alpha=0.9)
#ax3.plot(t_range, -envelope_imag_hilbert, color=pastel_orange, linewidth=2.5, alpha=0.5, linestyle='--')
for zt in zeros_t:
    if zt <= 100:
        ax3.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax3.axhline(y=0, color='black', linewidth=0.5)
ax3.set_xlabel('t', fontsize=12)
ax3.set_ylabel('Imaginary part', fontsize=12)
ax3.set_title('Imaginary Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 100)
ax3.set_ylim(-5, 5)

# Plot 4: Detail around selected zero
ax4 = plt.subplot(3, 3, 4)
detail_range = (t_range >= zeros_t[myzero] - 5) & (t_range <= zeros_t[myzero] + 5)
ax4.plot(t_range[detail_range], abs_standard[detail_range], color=pastel_blue, 
         linestyle='-', linewidth=2, label='|ζ(0.5+it)|', alpha=0.8)
ax4.plot(t_range[detail_range], abs_modified[detail_range], color=pastel_coral, 
         linestyle='--', linewidth=2, label='|ζ_mod(0.5+it)|', alpha=0.7)
ax4.plot(t_range[detail_range], envelope_abs_hilbert[detail_range], color=pastel_orange, 
         linewidth=3, label='Envelope curve', alpha=0.9)
ax4.axvline(x=zeros_t[myzero], color='green', linestyle='--', linewidth=2, label=f'Zero #{myzero+1}')
ax4.set_xlabel('t', fontsize=12)
ax4.set_ylabel('Modulus', fontsize=12)
ax4.set_title(f'Detail around zero (t ≈ {zeros_t[myzero]:.2f})', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

# Plot 5: Real part - detail
ax5 = plt.subplot(3, 3, 5)
ax5.plot(t_range[detail_range], real_standard[detail_range], color=pastel_blue, 
         linestyle='-', linewidth=2, label='Re(ζ)', alpha=0.8)
ax5.plot(t_range[detail_range], real_modified[detail_range], color=pastel_coral, 
         linestyle='--', linewidth=2, label='Re(ζ_mod)', alpha=0.7)
#ax5.plot(t_range[detail_range], envelope_real_hilbert[detail_range], color=pastel_orange, 
#         linewidth=3, label='Envelope curve', alpha=0.9)
#ax5.plot(t_range[detail_range], -envelope_real_hilbert[detail_range], color=pastel_orange, 
#         linewidth=3, alpha=0.5, linestyle='--')
ax5.axvline(x=zeros_t[myzero], color='green', linestyle='--', linewidth=2)
ax5.axhline(y=0, color='black', linewidth=0.5)
ax5.set_xlabel('t', fontsize=12)
ax5.set_ylabel('Real part', fontsize=12)
ax5.set_title('Re - Detail around zero', fontsize=14, fontweight='bold')
ax5.legend(fontsize=10)
ax5.grid(True, alpha=0.3)

# Plot 6: Imaginary part - detail
ax6 = plt.subplot(3, 3, 6)
ax6.plot(t_range[detail_range], imag_standard[detail_range], color=pastel_blue, 
         linestyle='-', linewidth=2, label='Im(ζ)', alpha=0.8)
ax6.plot(t_range[detail_range], imag_modified[detail_range], color=pastel_coral, 
         linestyle='--', linewidth=2, label='Im(ζ_mod)', alpha=0.7)
#ax6.plot(t_range[detail_range], envelope_imag_hilbert[detail_range], color=pastel_orange, 
#         linewidth=3, label='Envelope curve', alpha=0.9)
#ax6.plot(t_range[detail_range], -envelope_imag_hilbert[detail_range], color=pastel_orange, 
#         linewidth=3, alpha=0.5, linestyle='--')
ax6.axvline(x=zeros_t[myzero], color='green', linestyle='--', linewidth=2)
ax6.axhline(y=0, color='black', linewidth=0.5)
ax6.set_xlabel('t', fontsize=12)
ax6.set_ylabel('Imaginary part', fontsize=12)
ax6.set_title('Im - Detail around zero', fontsize=14, fontweight='bold')
ax6.legend(fontsize=10)
ax6.grid(True, alpha=0.3)

# Plot 7: Absolute difference in modulus - SAFE LOG SCALE
ax7 = plt.subplot(3, 3, 7)
ax7.plot(t_range, diff_abs, color=pastel_purple, linewidth=1.5)
for zt in zeros_t:
    if zt <= 100:
        ax7.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax7.set_xlabel('t', fontsize=12)
ax7.set_ylabel('||ζ_mod| - |ζ||', fontsize=12)
ax7.set_title('Difference in Modulus', fontsize=14, fontweight='bold')
ax7.grid(True, alpha=0.3)
ax7.set_xlim(0, 100)
try:
    ax7.set_yscale('log')
except:
    print("Warning: Cannot set log scale for plot 7")

# Plot 8: Difference in real part - SAFE LOG SCALE  
ax8 = plt.subplot(3, 3, 8)
ax8.plot(t_range, diff_real, color=pastel_orange, linewidth=1.5)
for zt in zeros_t:
    if zt <= 100:
        ax8.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax8.set_xlabel('t', fontsize=12)
ax8.set_ylabel('|Re(ζ_mod) - Re(ζ)|', fontsize=12)
ax8.set_title('Difference in Real Part', fontsize=14, fontweight='bold')
ax8.grid(True, alpha=0.3)
ax8.set_xlim(0, 100)
try:
    ax8.set_yscale('log')
except:
    print("Warning: Cannot set log scale for plot 8")

# Plot 9: Complex distance - SAFE LOG SCALE
ax9 = plt.subplot(3, 3, 9)
ax9.plot(t_range, diff_complex, color=pastel_coral, linewidth=1.5)
for zt in zeros_t:
    if zt <= 100:
        ax9.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=0.3)
ax9.set_xlabel('t', fontsize=12)
ax9.set_ylabel('|ζ_mod(0.5+it) - ζ(0.5+it)|', fontsize=12)
ax9.set_title('Complex Distance (total error)', fontsize=14, fontweight='bold')
ax9.grid(True, alpha=0.3)
ax9.set_xlim(0, 100)
try:
    ax9.set_yscale('log')
except:
    print("Warning: Cannot set log scale for plot 9")

try:
    plt.tight_layout()
except:
    print("Warning: tight_layout failed, continuing without it")
    
#plt.savefig('zeta_critical_final.png', dpi=200, bbox_inches='tight')
plt.show()

Computing standard zeta function...
  Progress: 0/2000
  Progress: 200/2000
  Progress: 400/2000
  Progress: 600/2000
  Progress: 800/2000
  Progress: 1000/2000
  Progress: 1200/2000
  Progress: 1400/2000
  Progress: 1600/2000
  Progress: 1800/2000
Computing modified zeta function...
  Progress: 0/2000
  Progress: 200/2000
  Progress: 400/2000
  Progress: 600/2000
  Progress: 800/2000
  Progress: 1000/2000
  Progress: 1200/2000
  Progress: 1400/2000
  Progress: 1600/2000
  Progress: 1800/2000
Computing envelope curves...
  Finding local maxima for modulus...

Creating plots...

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

20. From the Eta Function to the Zeta Function: Riemann's Transformation and an Arithmetic Observation

Structure of the Two Versions

v2 (modified_zeta_v2) — doubled even eta:

$$\zeta_{v2}(s) = 2\pi^{s-1}\cdot 2^s \sum_{\substack{n=2 \\ n\,\text{even}}}^{N} \left(1 - n^{-s}\right)(-1)^{n/2}$$

Red curve. Contains exclusively even numbers $n = 2, 4, 6, \ldots$
The result is $2\pi^{s-1}\eta(s)$ plus an additive shift — a doubled version of the eta function's structure.

v3 (modified_zeta_v3) — normalisation via Riemann's transformation:

$$\zeta_{v3}(s) = \frac{\zeta_{v2}(s)}{1 - 2^{1-s}}$$

Purple curve. Riemann used the identity $\eta(s) = (1-2^{1-s})\zeta(s)$ to extend the zeta function to the strip $0 < \text{Re}(s) < 1$.

Key Observation: Division Does Not Introduce Odd Numbers

The factor $(1 - 2^{1-s})$ in the denominator is a purely algebraic operation on the sum already computed.
It adds no new terms to the series — odd numbers $n = 1, 3, 5, \ldots$ are absent both before and after the division:

$$\zeta_{v3}(s) = \underbrace{\frac{1.25 \cdot 2}{2^s}\cdot\pi^{s-1}}_{\text{scaling factor, no zeros}} \cdot \zeta(s)$$

Algebraically this is $\zeta(s)$ multiplied by a real factor that has no zeros on the critical line.
The zeros of $\zeta_{v3}$ are identical to those of $\zeta(s)$ — and they originate exclusively from even numbers.

Odd numbers $n = 1, 3, 5, \ldots$ are strictly zero throughout this construction:
$\sin\!\left(\frac{\pi(1-n)}{2}\right) = 0$ for every odd $n$.
Riemann's normalisation does not change this fact.

Convergence of Minima to the Known Zeros

For finite $N$ the minima of $|\zeta_{v3,N}|$ are approximate — as $N$ grows they converge to the exact zeros of $\zeta$:

| $N$ | $|\zeta_{v3}(0.5 + 14.135i)|$ | Zero location error | |-----|------|------| | 100 | ≈ 0.030 | $\sim 10^{-2}$ | | 500 | ≈ 0.013 | $\sim 10^{-3}$ | | 2000 | ≈ 0.007 | $\sim 10^{-3}$ | | 5000 | ≈ 0.004 | $\sim 10^{-3}$ | | $\infty$ | $0$ | exact zero |

The rate of convergence is $\sim 1/\sqrt{N}$ — typical for oscillatory sums at the boundary of convergence.
Even so, the first known zeros are already localised to within hundredths for $N = 100$, and to within thousandths for $N = 2000$.

$$\lim_{N\to\infty} \zeta_{v3,N}(s) = \frac{1.25 \cdot 2}{2^s}\cdot\pi^{s-1}\cdot\zeta(s)$$

import numpy as np
import matplotlib.pyplot as plt
from mpmath import zeta as mpmath_zeta
import warnings
from scipy.signal import hilbert, find_peaks
from scipy.interpolate import interp1d
from scipy.special import zeta, gamma
from scipy.ndimage import gaussian_filter1d
from mpmath import zeta as mpmath_zeta, gamma as mpmath_gamma
from math import sqrt, floor
import math

warnings.filterwarnings('ignore')
myzero = 23
chartlim = 100
# High resolution for critical line
t_range = np.linspace(0, chartlim, 2000)

def is_prime(n):
  for i in range(2,int(math.sqrt(n))+1):
    if (n%i) == 0:
      return False
  return True

# Define pastel colors
pastel_coral = '#FF6B6B'
pastel_blue = '#6BB6FF'
pastel_purple = '#B57EDC'
pastel_orange = '#FFB84D'

def modified_zeta_v2(s, n_terms=5000):
    total = complex(0, 0)
    
    for n in range(1, n_terms + 1):
        if n%2 == 0:
            total_plus2 = ((n**s - 1) / n**s) * np.sin(np.pi * (1 - n) / 2) * np.pi**(s-1) * 2**s
            total += 2 * total_plus2

    total = total

    return total

def modified_zeta_v3(s, n_terms=2000):
    total = complex(0, 0)
    
    #s_complement = 1 - s
    #gamma_val = complex(mpmath_gamma(s_complement))

    for n in range(2, n_terms + 1):

        if n%2 == 0:
            total_plus2 = ((n**s - 1) / n**s) * np.sin(np.pi * (1 - n) / 2) * np.pi**(s-1) * 1.25 #for s=0.5 1.25 for s=0.2 1.45 for s=0.7 1.15
            total += 2 * total_plus2
        #if n%2 == 1:
        #    total_plus3 = ((n**s - 1) / n**s) * np.sin(np.pi * (0 - n) / 2) * np.pi**(s-1) * 2**s
        #    total += 2 * total_plus3

    total = total / (1 - (1/2**(s-1)))

    return total


def riemann_zeta(s):
    """Standard Riemann zeta function"""
    try:
        return complex(mpmath_zeta(s))
    except:
        return np.nan + 1j*np.nan

print("Computing standard zeta function...")
zeta_standard = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = riemann_zeta(s)
    zeta_standard.append(z)

print("Computing modified zeta function...")
zeta_modified = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = modified_zeta_v2(s, n_terms=2000)
    zeta_modified.append(z)

zeta_modified3 = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = modified_zeta_v3(s, n_terms=2000)
    zeta_modified3.append(z)

# Extract components - as numpy arrays
real_standard = np.array([z.real for z in zeta_standard])
imag_standard = np.array([z.imag for z in zeta_standard])
abs_standard = np.array([abs(z) for z in zeta_standard])

real_modified = np.array([z.real for z in zeta_modified])
imag_modified = np.array([z.imag for z in zeta_modified])
abs_modified = np.array([abs(z) for z in zeta_modified])

real_modified3 = np.array([z.real for z in zeta_modified3])
imag_modified3 = np.array([z.imag for z in zeta_modified3])
abs_modified3 = np.array([abs(z) for z in zeta_modified3])
#print(*abs_modified3)
# Compute envelope curves
print("Computing envelope curves...")

# For real and imaginary parts: Hilbert transform
envelope_real_hilbert = np.abs(hilbert(real_modified))
envelope_imag_hilbert = np.abs(hilbert(imag_modified))

# For modulus: Local maxima + interpolation
print("  Finding local maxima for modulus...")
peaks_abs, _ = find_peaks(abs_modified, distance=10, prominence=0.1)

if len(peaks_abs) > 5:
    envelope_abs_interp = interp1d(
        t_range[peaks_abs], 
        abs_modified[peaks_abs],
        kind='cubic',
        fill_value='extrapolate'
    )
    envelope_abs_hilbert = envelope_abs_interp(t_range)
    envelope_abs_hilbert = gaussian_filter1d(envelope_abs_hilbert, sigma=5)
else:
    def moving_maximum(signal, window_size):
        envelope = np.zeros_like(signal)
        half_window = window_size // 2
        for i in range(len(signal)):
            start = max(0, i - half_window)
            end = min(len(signal), i + half_window + 1)
            envelope[i] = np.max(signal[start:end])
        return envelope
    envelope_abs_hilbert = moving_maximum(abs_modified, window_size=50)

# Differences - IMPORTANT: Add epsilon and remove NaN/inf!
epsilon = 1e-10
diff_real_raw = np.abs(real_modified - real_standard)
diff_imag_raw = np.abs(imag_modified - imag_standard)
diff_abs_raw = np.abs(abs_modified - abs_standard)
diff_complex_raw = np.array([abs(zeta_modified[i] - zeta_standard[i]) for i in range(len(t_range))])

# Replace NaN and inf with epsilon values
diff_real = np.where(np.isfinite(diff_real_raw), diff_real_raw, epsilon) + epsilon
diff_imag = np.where(np.isfinite(diff_imag_raw), diff_imag_raw, epsilon) + epsilon
diff_abs = np.where(np.isfinite(diff_abs_raw), diff_abs_raw, epsilon) + epsilon
diff_complex = np.where(np.isfinite(diff_complex_raw), diff_complex_raw, epsilon) + epsilon

# First nontrivial zeros
zeros_t = [14.134725, 21.022040, 25.010858, 30.424876, 32.935062, 37.586178,
           40.918719, 43.327073, 48.005151, 49.773832, 52.970321, 56.446248,
           59.347044, 60.831779, 65.112544, 67.079810, 69.546402, 72.067158,
           75.704691, 77.144840, 79.337375, 82.910381, 84.735493, 87.425275,
           88.809111, 92.491899, 94.651344, 95.870634, 98.831194]

xticks = range(0, int(chartlim) + 1, 50)

fig2, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, abs_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='|ζ(0.5+it)|', alpha=1.0)
ax.plot(t_range, abs_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='|ζ_mod(0.5+it)|', alpha=1.0)
ax.plot(t_range, abs_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='|ζ_mod3(s+it)|', alpha=1.0)
#ax.plot(t_range, envelope_abs_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve (Hilbert)', alpha=0.9)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Modulus', fontsize=12)
ax.set_title('Modulus on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(0, 5)
#plt.savefig('zeta_modified.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()


fig3, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, real_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Re(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, real_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Re(ζ_mod(0.5+it))', alpha=1.0)
ax.plot(t_range, real_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='Re(ζ_mod3(s+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Real part', fontsize=12)
ax.set_title('Real Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedr.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

fig4, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, imag_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Im(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, imag_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Im(ζ_mod(0.5+it))', alpha=1.0)
ax.plot(t_range, imag_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='Im(ζ_mod3(s+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Imaginary part', fontsize=12)
ax.set_title('Imaginary Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedi.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

Computing standard zeta function...
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
Computing modified zeta function...
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
Computing envelope curves...
  Finding local maxima for modulus...

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

The Exact Coefficient $C(s)$ and the Definition of the $\delta$-Function

The coefficient in the code (1.25 for $s=0.5$, 1.45 for $s=0.2$, 1.15 for $s=0.7$) is not an empirical constant —
it is an exact function:

$$C(s) = \left(\frac{\pi}{2}\right)^{1-s}$$

Verification:

$s$	$C(s) = (\pi/2)^{1-s}$	value in code
0.2	$(\pi/2)^{0.8}$	≈ 1.4351 ≈ 1.45
0.5	$(\pi/2)^{0.5}$	≈ 1.2533 ≈ 1.25
0.7	$(\pi/2)^{0.3}$	≈ 1.1451 ≈ 1.15

The $\delta$-Function: First Definition

We introduce a new function $\delta(s)$ that encodes this choice of coefficient analytically:

$$\boxed{\delta(s) = \sum_{\substack{n=2 \\ n\,\text{even}}}^{\infty} 2\cdot\frac{n^s-1}{n^s}\cdot\sin\frac{\pi(1-n)}{2}\cdot\left(\frac{\pi}{2}\right)^{1-s}\cdot\frac{1}{1-2^{1-s}}}$$

Proof: $\delta(s) = \zeta(s)$

Expanding the factors algebraically:

$$\delta(s) = 2\cdot\underbrace{\left(\frac{\pi}{2}\right)^{1-s}}_{C(s)}\cdot\frac{\pi^{s-1}}{2^s}\cdot\frac{\eta(s)}{1-2^{1-s}}$$

$$= 2\cdot\frac{\pi^{1-s}}{2^{1-s}}\cdot\frac{\pi^{s-1}}{2^s}\cdot\zeta(s)$$

$$= 2\cdot\frac{\pi^{(1-s)+(s-1)}}{2^{(1-s)+s}}\cdot\zeta(s)$$

$$= 2\cdot\frac{\pi^{0}}{2^{1}}\cdot\zeta(s) = \zeta(s)$$

The powers of $\pi$ and $2$ cancel exactly — $\delta(s) \equiv \zeta(s)$ analytically, for all $s \in \mathbb{C}\setminus\{1\}$.
The code uses rounded values of $C(s)$; with the exact choice $(\pi/2)^{1-s}$, v3 reproduces the Riemann zeta function precisely.

How $C(s)$ Was Found

The coefficient $C(s)$ was not estimated numerically — it was derived analytically by requiring
that the entire expression $\delta(s)$ reduce exactly to $\zeta(s)$.

The condition $\delta(s) = \zeta(s)$ gives:

$$2 \cdot C(s) \cdot \frac{\pi^{s-1}}{2^s} \stackrel{!}{=} 1$$

$$C(s) = \frac{2^s}{2\cdot\pi^{s-1}} = \frac{2^{s-1}}{\pi^{s-1}} = \left(\frac{2}{\pi}\right)^{s-1} = \left(\frac{\pi}{2}\right)^{1-s}$$

The numerical values in the earlier code (1.25, 1.45, 1.15) were approximations of this expression
for specific values of $s$ — they retrospectively confirm the derivation.

Corrected implementation of modified_zeta_v3:

def modified_zeta_v3(s, n_terms=2000):
    """
    δ-function: exact representation of ζ(s) via even numbers
    δ(s) = Σ_{n even} 2·(n^s-1)/n^s · sin(π(1-n)/2) · (π/2)^(1-s) · 1/(1-2^(1-s))
    Algebraically: δ(s) = ζ(s)
    """
    total = complex(0, 0)
    C = (np.pi / 2) ** (1 - s.real)     # exact coefficient C(s) = (π/2)^(1-s)

    for n in range(2, n_terms + 1):
        if n % 2 == 0:
            total += 2 * ((n**s - 1) / n**s) * np.sin(np.pi * (1 - n) / 2) * np.pi**(s-1) * C

    total = total / (1 - 2**(1 - s))    # Riemann normalisation: η → ζ

    return total

Validity of $\delta(s) = \zeta(s)$ for the Entire Strip $\text{Re}(s) \in (0, 1)$

The algebraic proof does not require $\sigma = \frac{1}{2}$ — the identity holds for all
$s = \sigma + it$ with $\sigma \in (0,1)$, wherever $\eta(s)$ converges, i.e. for $\text{Re}(s) > 0$.

The key is that $C(s) = (\pi/2)^{1-s}$ depends only on $\sigma = \text{Re}(s)$
and the powers of $\pi$ and $2$ cancel for arbitrary $s$:

$$\delta(s) = 2\cdot\left(\frac{\pi}{2}\right)^{1-s}\cdot\frac{\pi^{s-1}}{2^s}\cdot\zeta(s) = \frac{2\pi^0}{2^1}\cdot\zeta(s) = \zeta(s) \quad \forall\, \text{Re}(s)\in(0,1)$$

For real $s \in (0,1)$ numerical convergence for finite $N$ is excellent
and improves monotonically with both $\sigma$ and $N$:

$s$ (real)	$\delta(s)$, $N=2000$	exact $\zeta(s)$	rel. error
0.5	−1.42219	−1.46035	2.6 %
0.7	−2.76121	−2.77839	0.6 %
0.8	−4.42416	−4.43754	0.3 %
0.9	−9.41622	−9.43011	0.1 %

For $s=0.8$ the error decreases as $\sim 1/N$:
$N=100$: 3.3 %, $N=1000$: 0.5 %, $N=10000$: 0.08 %.

Note on numerical convergence: For complex $s = \sigma + it$ with large $|t|$,
phase oscillations $e^{-it\ln n}$ dominate and slow the convergence of the finite sum
to its analytic limit. The approximation of zero locations on the critical line works well even for small $N$
due to oscillatory cancellation — but this is a property specific to $\sigma = \frac{1}{2}$.
The analytic identity $\delta(s) = \zeta(s)$ holds exactly in the limit $N \to \infty$ throughout the strip.

import numpy as np
import matplotlib.pyplot as plt
from mpmath import zeta as mpmath_zeta
import warnings
from scipy.signal import hilbert, find_peaks
from scipy.interpolate import interp1d
from scipy.special import zeta, gamma
from scipy.ndimage import gaussian_filter1d
from mpmath import zeta as mpmath_zeta, gamma as mpmath_gamma
from math import sqrt, floor
import math

warnings.filterwarnings('ignore')
myzero = 23
chartlim = 100
# High resolution for critical line
t_range = np.linspace(0, chartlim, 2000)

def is_prime(n):
  for i in range(2,int(math.sqrt(n))+1):
    if (n%i) == 0:
      return False
  return True

# Define pastel colors
pastel_coral = '#FF6B6B'
pastel_blue = '#6BB6FF'
pastel_purple = '#B57EDC'
pastel_orange = '#FFB84D'

def modified_zeta_v2(s, n_terms=5000):
    total = complex(0, 0)
    
    for n in range(1, n_terms + 1):
        if n%2 == 0:
            total_plus2 = ((n**s - 1) / n**s) * np.sin(np.pi * (1 - n) / 2) * np.pi**(s-1) * 2**s
            total += 2 * total_plus2

    total = total

    return total

def modified_zeta_v3(s, n_terms=2000):
    """
    δ-funkce: přesná reprezentace ζ(s) přes sudá čísla
    δ(s) = Σ_{n sudé} 2·(n^s-1)/n^s · sin(π(1-n)/2) · (π/2)^(1-s) · 1/(1-2^(1-s))
    Algebraicky: δ(s) = ζ(s)
    """
    total = complex(0, 0)
    C = (np.pi / 2) ** (1 - s.real)     # přesný koeficient C(s) = (π/2)^(1-s)

    for n in range(2, n_terms + 1):
        if n % 2 == 0:
            total += 2 * ((n**s - 1) / n**s) * np.sin(np.pi * (1 - n) / 2) * np.pi**(s-1) * C

    total = total / (1 - 2**(1 - s))    # Riemannova normalizace: η → ζ

    return total


def riemann_zeta(s):
    """Standard Riemann zeta function"""
    try:
        return complex(mpmath_zeta(s))
    except:
        return np.nan + 1j*np.nan

print("Computing standard zeta function...")
zeta_standard = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.8, t)
    z = riemann_zeta(s)
    zeta_standard.append(z)

print("Computing modified zeta function...")
zeta_modified = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.8, t)
    z = modified_zeta_v2(s, n_terms=2000)
    zeta_modified.append(z)

zeta_modified3 = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.8, t)
    z = modified_zeta_v3(s, n_terms=2000)
    zeta_modified3.append(z)

# Extract components - as numpy arrays
real_standard = np.array([z.real for z in zeta_standard])
imag_standard = np.array([z.imag for z in zeta_standard])
abs_standard = np.array([abs(z) for z in zeta_standard])

real_modified = np.array([z.real for z in zeta_modified])
imag_modified = np.array([z.imag for z in zeta_modified])
abs_modified = np.array([abs(z) for z in zeta_modified])

real_modified3 = np.array([z.real for z in zeta_modified3])
imag_modified3 = np.array([z.imag for z in zeta_modified3])
abs_modified3 = np.array([abs(z) for z in zeta_modified3])
#print(*abs_modified3)
# Compute envelope curves
print("Computing envelope curves...")

# For real and imaginary parts: Hilbert transform
envelope_real_hilbert = np.abs(hilbert(real_modified))
envelope_imag_hilbert = np.abs(hilbert(imag_modified))

# For modulus: Local maxima + interpolation
print("  Finding local maxima for modulus...")
peaks_abs, _ = find_peaks(abs_modified, distance=10, prominence=0.1)

if len(peaks_abs) > 5:
    envelope_abs_interp = interp1d(
        t_range[peaks_abs], 
        abs_modified[peaks_abs],
        kind='cubic',
        fill_value='extrapolate'
    )
    envelope_abs_hilbert = envelope_abs_interp(t_range)
    envelope_abs_hilbert = gaussian_filter1d(envelope_abs_hilbert, sigma=5)
else:
    def moving_maximum(signal, window_size):
        envelope = np.zeros_like(signal)
        half_window = window_size // 2
        for i in range(len(signal)):
            start = max(0, i - half_window)
            end = min(len(signal), i + half_window + 1)
            envelope[i] = np.max(signal[start:end])
        return envelope
    envelope_abs_hilbert = moving_maximum(abs_modified, window_size=50)

# Differences - IMPORTANT: Add epsilon and remove NaN/inf!
epsilon = 1e-10
diff_real_raw = np.abs(real_modified - real_standard)
diff_imag_raw = np.abs(imag_modified - imag_standard)
diff_abs_raw = np.abs(abs_modified - abs_standard)
diff_complex_raw = np.array([abs(zeta_modified[i] - zeta_standard[i]) for i in range(len(t_range))])

# Replace NaN and inf with epsilon values
diff_real = np.where(np.isfinite(diff_real_raw), diff_real_raw, epsilon) + epsilon
diff_imag = np.where(np.isfinite(diff_imag_raw), diff_imag_raw, epsilon) + epsilon
diff_abs = np.where(np.isfinite(diff_abs_raw), diff_abs_raw, epsilon) + epsilon
diff_complex = np.where(np.isfinite(diff_complex_raw), diff_complex_raw, epsilon) + epsilon

# First nontrivial zeros
zeros_t = [14.134725, 21.022040, 25.010858, 30.424876, 32.935062, 37.586178,
           40.918719, 43.327073, 48.005151, 49.773832, 52.970321, 56.446248,
           59.347044, 60.831779, 65.112544, 67.079810, 69.546402, 72.067158,
           75.704691, 77.144840, 79.337375, 82.910381, 84.735493, 87.425275,
           88.809111, 92.491899, 94.651344, 95.870634, 98.831194]

xticks = range(0, int(chartlim) + 1, 50)

fig2, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, abs_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='|ζ(0.8+it)|', alpha=1.0)
ax.plot(t_range, abs_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='|ζ_mod(0.8+it)|', alpha=1.0)
ax.plot(t_range, abs_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='|ζ_mod3(0.8+it)|', alpha=1.0)
#ax.plot(t_range, envelope_abs_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve (Hilbert)', alpha=0.9)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Modulus', fontsize=12)
ax.set_title('Modulus on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(0, 5)
#plt.savefig('zeta_modified.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()


fig3, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, real_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Re(ζ(0.8+it))', alpha=1.0)
ax.plot(t_range, real_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Re(ζ_mod(0.8+it))', alpha=1.0)
ax.plot(t_range, real_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='Re(ζ_mod3(0.8+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Real part', fontsize=12)
ax.set_title('Real Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedr.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

fig4, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, imag_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Im(ζ(0.8+it))', alpha=1.0)
ax.plot(t_range, imag_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Im(ζ_mod(0.8+it))', alpha=1.0)
ax.plot(t_range, imag_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='Im(ζ_mod3(0.8+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Imaginary part', fontsize=12)
ax.set_title('Imaginary Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedi.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

Computing standard zeta function...
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
Computing modified zeta function...
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
Computing envelope curves...
  Finding local maxima for modulus...

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

21. The Dirichlet Beta Function: Moving from Even to Odd Numbers

Shifting the Sine by 1: Selecting Odd Numbers

Replacing the sine argument from $\sin\!\left(\frac{\pi(1-n)}{2}\right)$ to $\sin\!\left(\frac{\pi(0-n)}{2}\right)$
switches the sieve — even terms are now zero and odd $n$ are retained:

$n$	$\sin\!\left(\frac{\pi(1-n)}{2}\right)$ — even (zeta)	$\sin\!\left(\frac{\pi(0-n)}{2}\right)$ — odd (beta)
1	0	−1
2	−1	0
3	0	+1
4	+1	0
5	0	−1
6	−1	0
7	0	+1
8	+1	0

For odd $n = 2k-1$ we have $\sin\!\left(\frac{-\pi n}{2}\right) = (-1)^k$,
so the signs run $-1, +1, -1, +1, \ldots$ for $n = 1, 3, 5, 7, \ldots$

This is a shift of the argument by 1 — otherwise the structure is identical to the $\delta$-function for even numbers.

The Dirichlet Beta Function

The alternating series over odd numbers defines the Dirichlet beta function:

$$\beta(s) = \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)^s} = \frac{1}{1^s} - \frac{1}{3^s} + \frac{1}{5^s} - \frac{1}{7^s} + \cdots$$

Just as the $\delta$-function reproduces $\zeta(s)$ via even numbers,
the red curve here approximates $\beta(s)$ via odd numbers.

Key Relationship: $\zeta \leftrightarrow \pi$, $\beta \leftrightarrow G$

Both functions are deeply connected to fundamental constants:

Function	Special values	Constant
$\zeta(2) = \frac{\pi^2}{6}$	$\zeta(4) = \frac{\pi^4}{90}$, $\zeta(6) = \frac{\pi^6}{945}$, $\ldots$	$\pi$
$\beta(1) = \frac{\pi}{4}$	$\beta(3) = \frac{\pi^3}{32}$, $\beta(5) = \frac{5\pi^5}{1536}$, $\ldots$	$\pi$
$\beta(2) = G$	$\beta(4) = \frac{\pi^4 - ?}{?}$ (open problem)	Catalan's constant

The zeta function at even arguments always yields rational multiples of powers of $\pi$.
The beta function at odd arguments does likewise — but $\beta(2) = G$ stands outside this pattern:
Catalan's constant $G \approx 0.9159655942$ has so far not been expressed in terms of $\pi$,
nor has it been proved irrational.

This gap — $G$ as a value of $\beta$ at an even argument with no closed form —
is the subject of the next cell, where we bring both functions together.

The Purple Curve: Products of Odd Sums and Prime Numbers

modified_zeta_v3 in this block computes $(\text{odd sum})^2$ — the square of the sum over odd numbers.
The result is the characteristic arc shape, because multiplying a complex number by itself
doubles the phase and squares the modulus.

A hypothetical route to prime numbers: numbers with no non-trivial divisors (primes) should produce
sharp edges in such a sum, while composite numbers (multiples of earlier primes)
would yield rounded arcs reflecting their factor structure.
This direction requires replacing the squaring $(\cdot)^2$ with a summation function over divisors.

import numpy as np
import matplotlib.pyplot as plt
from mpmath import zeta as mpmath_zeta
import warnings
from scipy.signal import hilbert, find_peaks
from scipy.interpolate import interp1d
from scipy.special import zeta, gamma
from scipy.ndimage import gaussian_filter1d
from mpmath import zeta as mpmath_zeta, gamma as mpmath_gamma
from math import sqrt, floor
import math

warnings.filterwarnings('ignore')
myzero = 23
chartlim = 100
# High resolution for critical line
t_range = np.linspace(0, chartlim, 2000)

def is_prime(n):
  for i in range(2,int(math.sqrt(n))+1):
    if (n%i) == 0:
      return False
  return True

# Define pastel colors
pastel_coral = '#FF6B6B'
pastel_blue = '#6BB6FF'
pastel_purple = '#B57EDC'
pastel_orange = '#FFB84D'

def modified_zeta_v2(s, n_terms=2000):
    """
    δ-funkce: přesná reprezentace ζ(s) přes lichá čísla
    δ(s) = Σ_{n liché} 2·(n^s-1)/n^s · sin(π(0-n)/2) · (π/2)^(1-s) · 1/(1-2^(1-s))
    Algebraicky: δ(s) = ζ(s)
    """
    total = complex(0, 0)
    C = (np.pi / 2) ** (1 - s.real)     # přesný koeficient C(s) = (π/2)^(1-s)

    for n in range(2, n_terms + 1):
        if n % 2 == 1:
            total += 2 * ((n**s - 1) / n**s) * np.sin(np.pi * (0 - n) / 2) * np.pi**(s-1) * C

    total = total

    return total

def modified_zeta_v3(s, n_terms=2000):
    total = complex(0, 0)
    total2 = complex(0, 0)

    for n in range(2, n_terms + 1):

        if n%2 == 0:
            total_plus = ((n**s - 1) / n**s) * np.sin(np.pi * (1 - n) / 2) * np.pi**(s-1) * 1.25 #for s=0.5 1.25 for s=0.2 1.45 for s=0.7 1.15
            total += total_plus
        if n%2 == 1:
            total_plus2 = ((n**s - 1) / n**s) * np.sin(np.pi * (0 - n) / 2) * np.pi**(s-1) * 1.25
            total2 += total_plus2

    total = total2 * total2

    return total


def riemann_zeta(s):
    """Standard Riemann zeta function"""
    try:
        return complex(mpmath_zeta(s))
    except:
        return np.nan + 1j*np.nan

print("Computing standard zeta function...")
zeta_standard = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = riemann_zeta(s)
    zeta_standard.append(z)

print("Computing modified zeta function...")
zeta_modified = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = modified_zeta_v2(s, n_terms=2000)
    zeta_modified.append(z)

zeta_modified3 = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = modified_zeta_v3(s, n_terms=2000)
    zeta_modified3.append(z)

# Extract components - as numpy arrays
real_standard = np.array([z.real for z in zeta_standard])
imag_standard = np.array([z.imag for z in zeta_standard])
abs_standard = np.array([abs(z) for z in zeta_standard])

real_modified = np.array([z.real for z in zeta_modified])
imag_modified = np.array([z.imag for z in zeta_modified])
abs_modified = np.array([abs(z) for z in zeta_modified])

real_modified3 = np.array([z.real for z in zeta_modified3])
imag_modified3 = np.array([z.imag for z in zeta_modified3])
abs_modified3 = np.array([abs(z) for z in zeta_modified3])
#print(*abs_modified3)
# Compute envelope curves
print("Computing envelope curves...")

# For real and imaginary parts: Hilbert transform
envelope_real_hilbert = np.abs(hilbert(real_modified))
envelope_imag_hilbert = np.abs(hilbert(imag_modified))

# For modulus: Local maxima + interpolation
print("  Finding local maxima for modulus...")
peaks_abs, _ = find_peaks(abs_modified, distance=10, prominence=0.1)

if len(peaks_abs) > 5:
    envelope_abs_interp = interp1d(
        t_range[peaks_abs], 
        abs_modified[peaks_abs],
        kind='cubic',
        fill_value='extrapolate'
    )
    envelope_abs_hilbert = envelope_abs_interp(t_range)
    envelope_abs_hilbert = gaussian_filter1d(envelope_abs_hilbert, sigma=5)
else:
    def moving_maximum(signal, window_size):
        envelope = np.zeros_like(signal)
        half_window = window_size // 2
        for i in range(len(signal)):
            start = max(0, i - half_window)
            end = min(len(signal), i + half_window + 1)
            envelope[i] = np.max(signal[start:end])
        return envelope
    envelope_abs_hilbert = moving_maximum(abs_modified, window_size=50)

# Differences - IMPORTANT: Add epsilon and remove NaN/inf!
epsilon = 1e-10
diff_real_raw = np.abs(real_modified - real_standard)
diff_imag_raw = np.abs(imag_modified - imag_standard)
diff_abs_raw = np.abs(abs_modified - abs_standard)
diff_complex_raw = np.array([abs(zeta_modified[i] - zeta_standard[i]) for i in range(len(t_range))])

# Replace NaN and inf with epsilon values
diff_real = np.where(np.isfinite(diff_real_raw), diff_real_raw, epsilon) + epsilon
diff_imag = np.where(np.isfinite(diff_imag_raw), diff_imag_raw, epsilon) + epsilon
diff_abs = np.where(np.isfinite(diff_abs_raw), diff_abs_raw, epsilon) + epsilon
diff_complex = np.where(np.isfinite(diff_complex_raw), diff_complex_raw, epsilon) + epsilon

# First nontrivial zeros
zeros_t = [14.134725, 21.022040, 25.010858, 30.424876, 32.935062, 37.586178,
           40.918719, 43.327073, 48.005151, 49.773832, 52.970321, 56.446248,
           59.347044, 60.831779, 65.112544, 67.079810, 69.546402, 72.067158,
           75.704691, 77.144840, 79.337375, 82.910381, 84.735493, 87.425275,
           88.809111, 92.491899, 94.651344, 95.870634, 98.831194]

xticks = range(0, int(chartlim) + 1, 50)

fig2, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, abs_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='|ζ(0.5+it)|', alpha=1.0)
ax.plot(t_range, abs_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='|ζ_mod(0.5+it)|', alpha=1.0)
ax.plot(t_range, abs_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='|ζ_mod3(s+it)|', alpha=1.0)
#ax.plot(t_range, envelope_abs_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve (Hilbert)', alpha=0.9)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Modulus', fontsize=12)
ax.set_title('Modulus on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(0, 5)
#plt.savefig('zeta_modified.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()


fig3, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, real_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Re(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, real_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Re(ζ_mod(0.5+it))', alpha=1.0)
ax.plot(t_range, real_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='Re(ζ_mod3(s+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Real part', fontsize=12)
ax.set_title('Real Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedr.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

fig4, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, imag_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Im(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, imag_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Im(ζ_mod(0.5+it))', alpha=1.0)
ax.plot(t_range, imag_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='Im(ζ_mod3(s+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Imaginary part', fontsize=12)
ax.set_title('Imaginary Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedi.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

Computing standard zeta function...
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
Computing modified zeta function...
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
Computing envelope curves...
  Finding local maxima for modulus...

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

import numpy as np
import matplotlib.pyplot as plt
from mpmath import zeta as mpmath_zeta
import warnings
from scipy.signal import hilbert, find_peaks
from scipy.interpolate import interp1d
from scipy.special import zeta, gamma
from scipy.ndimage import gaussian_filter1d
from mpmath import zeta as mpmath_zeta, gamma as mpmath_gamma
from mpl_toolkits.mplot3d import Axes3D

warnings.filterwarnings('ignore')
myzero = 9
max_x = 100

def delta_even(s, n_terms=2000):
    """
    δ-funkce přes sudá čísla → ζ(s)
    δ(s) = Σ_{n sudé} 2·(n^s-1)/n^s · sin(π(1-n)/2) · (π/2)^(1-s) · 1/(1-2^(1-s))
    """
    total = complex(0, 0)
    C = (np.pi / 2) ** (1 - s.real)
    for n in range(2, n_terms + 1):
        if n % 2 == 0:
            total += 2 * ((n**s - 1) / n**s) * np.sin(np.pi * (1 - n) / 2) * np.pi**(s-1) * C
    total = total / (1 - 2**(1 - s))
    return total

def delta_odd(s, n_terms=2000):
    """
    δ-funkce přes lichá čísla → β(s) (Dirichletova beta)
    β(s) = Σ_{n liché} 2·(n^s-1)/n^s · sin(π(0-n)/2) · (π/2)^(1-s)
    """
    total = complex(0, 0)
    C = (np.pi / 2) ** (1 - s.real)
    for n in range(1, n_terms + 1):
        if n % 2 == 1:
            total += 2 * ((n**s - 1) / n**s) * np.sin(np.pi * (0 - n) / 2) * np.pi**(s-1) * C
    return total


# Výpočet
print("Computing ζ(s)...")
zeta_standard = []
for i, t in enumerate(t_range):
    if i % 200 == 0: print(f"  {i}/{len(t_range)}")
    zeta_standard.append(riemann_zeta(complex(0.5, t)))

print("Computing δ_even(s) → ζ(s)...")
zeta_even = []
for i, t in enumerate(t_range):
    if i % 200 == 0: print(f"  {i}/{len(t_range)}")
    zeta_even.append(delta_even(complex(0.5, t), n_terms=2000))

print("Computing δ_odd(s) → β(s)...")
zeta_odd = []
for i, t in enumerate(t_range):
    if i % 200 == 0: print(f"  {i}/{len(t_range)}")
    zeta_odd.append(delta_odd(complex(0.5, t), n_terms=2000))

real_standard  = np.array([z.real for z in zeta_standard])
imag_standard  = np.array([z.imag for z in zeta_standard])
real_even      = np.array([z.real for z in zeta_even])
imag_even      = np.array([z.imag for z in zeta_even])
real_odd       = np.array([z.real for z in zeta_odd])
imag_odd       = np.array([z.imag for z in zeta_odd])

# ── 3D pohled ──────────────────────────────────────────────────────────────
fig_3d = plt.figure(figsize=(20, 7))

datasets_3d = [
    (real_standard, imag_standard, pastel_blue,   'ζ(0.5+it) — přesná',       'Standard Zeta'),
    (real_even,     imag_even,     pastel_coral,   'δ_even(0.5+it) → ζ(s)',    'δ-funkce sudá → ζ'),
    (real_odd,      imag_odd,      pastel_purple,  'δ_odd(0.5+it) → β(s)',     'δ-funkce lichá → β'),
]

for k, (re, im, color, label, title) in enumerate(datasets_3d):
    ax = fig_3d.add_subplot(1, 3, k+1, projection='3d')
    ax.plot(t_range, re, im, color=color, linewidth=1.5, label=label, alpha=0.8)
    for zt in zeros_t:
        if zt <= max_x:
            idx = np.argmin(np.abs(t_range - zt))
            ax.scatter([zt], [re[idx]], [im[idx]], color='green', s=50, alpha=0.6)
    ax.set_xlabel('t', fontsize=10, labelpad=8)
    ax.set_ylabel('Re', fontsize=10, labelpad=8)
    ax.set_zlabel('Im', fontsize=10, labelpad=8)
    ax.set_title(title, fontsize=12, fontweight='bold', pad=10)
    ax.legend(fontsize=9)
    ax.view_init(elev=20, azim=45)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
#plt.savefig('zeta_3d_rotations.png', dpi=200, bbox_inches='tight')
plt.show()

# ── 2D komplexní rovina ────────────────────────────────────────────────────
colors_2d = plt.cm.plasma(np.linspace(0, 1, len(t_range)))

datasets_2d = [
    (real_standard, imag_standard, 'ζ(s) — Standard Zeta'),
    (real_even,     imag_even,     'δ_even(s) → ζ(s) — sudá čísla'),
    (real_odd,      imag_odd,      'δ_odd(s) → β(s) — lichá čísla'),
]

for re, im, title in datasets_2d:
    fig, ax = plt.subplots(figsize=(10, 10))
    for i in range(len(t_range) - 1):
        ax.plot(re[i:i+2], im[i:i+2], color=colors_2d[i], linewidth=1.5, alpha=0.7)
    for zt in zeros_t:
        if zt <= max_x:
            idx = np.argmin(np.abs(t_range - zt))
            ax.scatter(re[idx], im[idx], color='red', s=100, marker='o',
                      edgecolors='black', linewidths=2, zorder=5)
            ax.annotate(f't≈{zt:.1f}', (re[idx], im[idx]),
                       fontsize=7, alpha=0.8, xytext=(5, 5),
                       textcoords='offset points')
    ax.axhline(y=0, color='black', linewidth=0.8, linestyle='--', alpha=0.5)
    ax.axvline(x=0, color='black', linewidth=0.8, linestyle='--', alpha=0.5)
    ax.set_xlabel('Re', fontsize=11)
    ax.set_ylabel('Im', fontsize=11)
    ax.set_title(title, fontsize=13, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')
    plt.tight_layout()
    plt.show()

Computing ζ(s)...
  0/2000
  200/2000
  400/2000
  600/2000
  800/2000
  1000/2000
  1200/2000
  1400/2000
  1600/2000
  1800/2000
Computing δ_even(s) → ζ(s)...
  0/2000
  200/2000
  400/2000
  600/2000
  800/2000
  1000/2000
  1200/2000
  1400/2000
  1600/2000
  1800/2000
Computing δ_odd(s) → β(s)...
  0/2000
  200/2000
  400/2000
  600/2000
  800/2000
  1000/2000
  1200/2000
  1400/2000
  1600/2000
  1800/2000

22. Unifying $\zeta$, $\eta$, $\beta$ with a Single $\delta$-Function

The Generalised $\delta$-Function with Parameter $k$

$$\boxed{ \delta(k,\, s,\, N) = \frac{1}{1 - 2^{1-s}} \sum_{\substack{n=2 \\ n\,\text{even}}}^{N} 2 \cdot \frac{n^s-1}{n^s} \cdot \sin\!\left(\frac{\pi(k-n)}{2}\right) \cdot \left(\frac{\pi}{2}\right)^{\!1-s} }$$

The parameter $k \in \mathbb{R}$ rotates the sine argument — and thereby switches between functions:

$k$	active $n$	$\delta(k, s, \infty)$
$k = 1$	even	$\zeta(s)$ — Riemann zeta
$k = 0$	odd	$\beta(s)$ — Dirichlet beta
$k = \frac{1}{2}$	all, weight $\frac{\sqrt{2}}{2}$	$\frac{\sqrt{2}}{2}\!\left(\beta(s) + \frac{\eta(s)}{2^s}\right)$

The Shift $k = \tfrac{1}{2}$: Sine Pattern and Activation of All $n$

For $k = \tfrac{1}{2}$, $\sin\!\left(\frac{\pi(0.5-n)}{2}\right)$ takes the value $\pm\frac{\sqrt{2}}{2}$ for every $n$:

$n$	1	2	3	4	5	6	7	8
value	$-\frac{\sqrt{2}}{2}$	$-\frac{\sqrt{2}}{2}$	$+\frac{\sqrt{2}}{2}$	$+\frac{\sqrt{2}}{2}$	$-\frac{\sqrt{2}}{2}$	$-\frac{\sqrt{2}}{2}$	$+\frac{\sqrt{2}}{2}$	$+\frac{\sqrt{2}}{2}$

Both odd and even numbers contribute with equal weight $\frac{\sqrt{2}}{2}$ — $\beta$ and $\eta$ are present symmetrically.
The sign pattern has period 4 instead of 2 — the transition $\theta_3 \to \theta_2$ in Jacobi theta functions.

The Rotation Identity

For arbitrary $k \in \mathbb{R}$ and $s \in \mathbb{C}$, in the limit $N \to \infty$:

$$\delta(k, s, \infty) = \beta(s)\cdot\cos\!\left(\frac{\pi k}{2}\right) + \frac{\eta(s)}{2^s}\cdot\sin\!\left(\frac{\pi k}{2}\right)$$

The shift $k$ rotates the vector $\bigl(\beta(s),\; \eta(s)/2^s\bigr)$ by the angle $\frac{\pi k}{2}$:

$$\begin{pmatrix} \delta(0,s) \\ \delta(1,s) \end{pmatrix} = \begin{pmatrix} \cos\frac{\pi k}{2} & -\sin\frac{\pi k}{2} \\[4pt] \sin\frac{\pi k}{2} & \cos\frac{\pi k}{2} \end{pmatrix} \begin{pmatrix} \beta(s) \\[4pt] \eta(s)/2^s \end{pmatrix}$$

The transition $k\!: 0 \to 1$ (swap) exchanges the roles of $\beta$ and $\eta/2^s$ — that is, it swaps odd and even numbers.

Special Values at $s = 2$: $G$ and $\pi^2/48$

$$\beta(2) = G = 0.9159655942\ldots \qquad \frac{\eta(2)}{2^2} = \frac{\pi^2}{48} = 0.2056167584\ldots \qquad \zeta(2) = \frac{\pi^2}{6}$$

The rotational swap at $s=2$:

$$\boxed{ \begin{pmatrix} \delta(0,2) \\ \delta(1,2) \end{pmatrix} = \underbrace{\begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}}_{M} \begin{pmatrix} G \\ \pi^2/48 \end{pmatrix} = \begin{pmatrix} -\pi^2/48 \\ G \end{pmatrix} }$$

$G$ and $\pi^2/48$ are components of a single 2D vector — the matrix $M$ has integer entries,
and rotation by $\frac{\pi}{2}$ is an automorphism of the same algebraic structure.

A Response to Grothendieck's Period Conjecture

Both functions — $\beta$ (odd $n$) and $\eta/2^s$ (even $n$) — carry identical weights
$\left(1 - n^{-s}\right)$, differing only in the rotation of the sine argument.

The Grothendieck–Kontsevich–Zagier conjecture asserts that algebraic relations
between periods arise solely from linearity, variable substitution, and Stokes's theorem.

Rotation by $\frac{\pi k}{2}$ is multiplication of the character $\chi_4$ by $i$
in the complex group of characters — an automorphism of the same structure,
not a geometrically independent operation. $G$ and $\pi^2$ are therefore not
algebraically independent under this rotation — they are components of a single vector.

Connection via the Polylogarithm $\mathrm{Li}_2$

Key identity — verified to 10 decimal places:

$$\mathrm{Li}_2(i) = -\frac{\pi^2}{48} + i\,G$$

The real part is $-\pi^2/48$ and the imaginary part is $G$ — precisely the components of the rotation vector.
Multiplication by $i$ in the complex plane corresponds to the swap $k\!: 0 \to 1$:

$$i \cdot \mathrm{Li}_2(i) = -G - i\,\frac{\pi^2}{48}$$

Abel's functional equation at the point $1+i$:

$$\mathrm{Im}\!\left(\mathrm{Li}_2(1+i)\right) = G + \frac{\pi \ln 2}{4} = 1.4603621168\ldots$$

This relation is a candidate for a proof of the transcendence of $G$ via Nesterenko's theorem,
which establishes the algebraic independence of $\pi$, $e^\pi$, and $\Gamma(1/4)$.

import numpy as np
import matplotlib.pyplot as plt
from mpmath import zeta as mpmath_zeta
import warnings
from scipy.signal import hilbert, find_peaks
from scipy.interpolate import interp1d
from scipy.special import zeta, gamma
from scipy.ndimage import gaussian_filter1d
from mpmath import zeta as mpmath_zeta, gamma as mpmath_gamma
from math import sqrt, floor
import math

warnings.filterwarnings('ignore')
myzero = 23
chartlim = 100
# High resolution for critical line
t_range = np.linspace(0, chartlim, 2000)


# Define pastel colors
pastel_coral = '#FF6B6B'
pastel_blue = '#6BB6FF'
pastel_purple = '#B57EDC'
pastel_orange = '#FFB84D'

def modified_zeta_v2(s, n_terms=2000):
    """
    δ-funkce: ζ(s) přes lichá čísla
    δ(s) = Σ_{n liché} 2·(n^s-1)/n^s · sin(π(0-n)/2) · (π/2)^(0.5-s) · 1/(1-2^(1-s))
    Algebraicky: δ(s) = ζ(s)
    """
    total = complex(0, 0)
    C = (np.pi / 2) ** (1 - s.real)     # přesný koeficient C(s) = (π/2)^(1-s)

    for n in range(2, n_terms + 1):
        if n % 2 == 1:
            total += 2 * ((n**s - 1) / n**s) * np.sin(np.pi * (0.5- n) / 2) * np.pi**(s-1) * C

    total = total

    return total

def modified_zeta_v3(s, n_terms=2000):
    """
    δ-funkce: ζ(s) přes sudá čísla
    δ(s) = Σ_{n sudé} 2·(n^s-1)/n^s · sin(π(0-n)/2) · (π/2)^(0.5-s) · 1/(1-2^(1-s))
    Algebraicky: δ(s) = ζ(s)
    """
    total = complex(0, 0)
    C = (np.pi / 2) ** (1 - s.real)     # přesný koeficient C(s) = (π/2)^(1-s)

    for n in range(2, n_terms + 1):
        if n % 2 == 0:
            total += 2 * ((n**s - 1) / n**s) * np.sin(np.pi * (0.5- n) / 2) * np.pi**(s-1) * C

    total = total

    return total


def riemann_zeta(s):
    """Standard Riemann zeta function"""
    try:
        return complex(mpmath_zeta(s))
    except:
        return np.nan + 1j*np.nan

print("Computing standard zeta function...")
zeta_standard = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = riemann_zeta(s)
    zeta_standard.append(z)

print("Computing modified zeta function...")
zeta_modified = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = modified_zeta_v2(s, n_terms=2000)
    zeta_modified.append(z)

zeta_modified3 = []
for i, t in enumerate(t_range):
    if i % chartlim == 0:
        print(f"  Progress: {i}/{len(t_range)}")
    s = complex(0.5, t)
    z = modified_zeta_v3(s, n_terms=2000)
    zeta_modified3.append(z)

# Extract components - as numpy arrays
real_standard = np.array([z.real for z in zeta_standard])
imag_standard = np.array([z.imag for z in zeta_standard])
abs_standard = np.array([abs(z) for z in zeta_standard])

real_modified = np.array([z.real for z in zeta_modified])
imag_modified = np.array([z.imag for z in zeta_modified])
abs_modified = np.array([abs(z) for z in zeta_modified])

real_modified3 = np.array([z.real for z in zeta_modified3])
imag_modified3 = np.array([z.imag for z in zeta_modified3])
abs_modified3 = np.array([abs(z) for z in zeta_modified3])
#print(*abs_modified3)
# Compute envelope curves
print("Computing envelope curves...")

# For real and imaginary parts: Hilbert transform
envelope_real_hilbert = np.abs(hilbert(real_modified))
envelope_imag_hilbert = np.abs(hilbert(imag_modified))

# For modulus: Local maxima + interpolation
print("  Finding local maxima for modulus...")
peaks_abs, _ = find_peaks(abs_modified, distance=10, prominence=0.1)

if len(peaks_abs) > 5:
    envelope_abs_interp = interp1d(
        t_range[peaks_abs], 
        abs_modified[peaks_abs],
        kind='cubic',
        fill_value='extrapolate'
    )
    envelope_abs_hilbert = envelope_abs_interp(t_range)
    envelope_abs_hilbert = gaussian_filter1d(envelope_abs_hilbert, sigma=5)
else:
    def moving_maximum(signal, window_size):
        envelope = np.zeros_like(signal)
        half_window = window_size // 2
        for i in range(len(signal)):
            start = max(0, i - half_window)
            end = min(len(signal), i + half_window + 1)
            envelope[i] = np.max(signal[start:end])
        return envelope
    envelope_abs_hilbert = moving_maximum(abs_modified, window_size=50)

# Differences - IMPORTANT: Add epsilon and remove NaN/inf!
epsilon = 1e-10
diff_real_raw = np.abs(real_modified - real_standard)
diff_imag_raw = np.abs(imag_modified - imag_standard)
diff_abs_raw = np.abs(abs_modified - abs_standard)
diff_complex_raw = np.array([abs(zeta_modified[i] - zeta_standard[i]) for i in range(len(t_range))])

# Replace NaN and inf with epsilon values
diff_real = np.where(np.isfinite(diff_real_raw), diff_real_raw, epsilon) + epsilon
diff_imag = np.where(np.isfinite(diff_imag_raw), diff_imag_raw, epsilon) + epsilon
diff_abs = np.where(np.isfinite(diff_abs_raw), diff_abs_raw, epsilon) + epsilon
diff_complex = np.where(np.isfinite(diff_complex_raw), diff_complex_raw, epsilon) + epsilon

# First nontrivial zeros
zeros_t = [14.134725, 21.022040, 25.010858, 30.424876, 32.935062, 37.586178,
           40.918719, 43.327073, 48.005151, 49.773832, 52.970321, 56.446248,
           59.347044, 60.831779, 65.112544, 67.079810, 69.546402, 72.067158,
           75.704691, 77.144840, 79.337375, 82.910381, 84.735493, 87.425275,
           88.809111, 92.491899, 94.651344, 95.870634, 98.831194]

xticks = range(0, int(chartlim) + 1, 50)

fig2, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, abs_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='|ζ(0.5+it)|', alpha=1.0)
ax.plot(t_range, abs_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='|ζ_mod(0.5+it)|', alpha=1.0)
ax.plot(t_range, abs_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='|ζ_mod3(s+it)|', alpha=1.0)
#ax.plot(t_range, envelope_abs_hilbert, color=pastel_orange, linewidth=2.5, label='Envelope curve (Hilbert)', alpha=0.9)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Modulus', fontsize=12)
ax.set_title('Modulus on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(0, 5)
#plt.savefig('zeta_modified.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()


fig3, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, real_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Re(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, real_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Re(ζ_mod(0.5+it))', alpha=1.0)
ax.plot(t_range, real_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='Re(ζ_mod3(s+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Real part', fontsize=12)
ax.set_title('Real Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedr.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

fig4, ax = plt.subplots(figsize=(chartlim, 6))
ax.plot(t_range, imag_standard, color=pastel_blue, linestyle='-', linewidth=1, 
        label='Im(ζ(0.5+it))', alpha=1.0)
ax.plot(t_range, imag_modified, color=pastel_coral, linestyle='--', linewidth=1, 
        label='Im(ζ_mod(0.5+it))', alpha=1.0)
ax.plot(t_range, imag_modified3, color=pastel_purple, linestyle='--', linewidth=1, 
        label='Im(ζ_mod3(s+it))', alpha=1.0)
for zt in zeros_t:
    if zt <= chartlim:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.5, alpha=1.0)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xticks(xticks)
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('Imaginary part', fontsize=12)
ax.set_title('Imaginary Part on Critical Line (overview)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, chartlim)
ax.set_ylim(-5, 5)
#plt.savefig('zeta_modifiedi.png', dpi=300, bbox_inches='tight')
plt.show()
plt.tight_layout()

Computing standard zeta function...
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
Computing modified zeta function...
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
  Progress: 0/2000
  Progress: 100/2000
  Progress: 200/2000
  Progress: 300/2000
  Progress: 400/2000
  Progress: 500/2000
  Progress: 600/2000
  Progress: 700/2000
  Progress: 800/2000
  Progress: 900/2000
  Progress: 1000/2000
  Progress: 1100/2000
  Progress: 1200/2000
  Progress: 1300/2000
  Progress: 1400/2000
  Progress: 1500/2000
  Progress: 1600/2000
  Progress: 1700/2000
  Progress: 1800/2000
  Progress: 1900/2000
Computing envelope curves...
  Finding local maxima for modulus...

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

<Figure size 640x480 with 0 Axes>

23. Rotation by the Parameter $k$: Continuously from $\beta$ to $\zeta$

Both functions in the previous block (modified_zeta_v2, modified_zeta_v3) are special cases
of a single generalised $\delta$-function with parameter $k \in [0,1]$:

$$\boxed{\delta(k, s, N) = \sum_{n=1}^{N} 2\cdot\frac{n^s-1}{n^s} \cdot\sin!\left(\frac{\pi(k-n)}{2}\right) \cdot\left(\frac{\pi}{2}\right)^{!1-s}}$$

The parameter $k$ rotates the sine argument and switches between three functions:

$k$	active $n$	sine weight	$\delta(k, s, \infty)$
$0$	odd	$\pm 1$	$\beta(s)$ — Catalan $G$
$\tfrac{1}{2}$	all	$\pm\tfrac{\sqrt{2}}{2}$	$\tfrac{\sqrt{2}}{2}\!\left(\beta + \tfrac{\eta}{2^s}\right)$
$1$	even	$\pm 1$	$\zeta(s)$ — $\pi^2$

Alternative Visualisation: Continuous Rotation in $(k, t)$ Space

Plot 1 shows three curves $|\delta(k, \frac{1}{2}+it)|$ for $k = 0, \frac{1}{2}, 1$
alongside the exact $\zeta$ (black). The zero locations (green vertical lines) migrate with $k$ —
for $k=0$ they are the zeros of $\beta$, for $k=1$ the zeros of $\zeta$.

Plot 2 (heatmap) displays $|\delta(k, \frac{1}{2}+it)|$ across the full $(k, t)$ space:

• $x$-axis: imaginary part $t \in [0, 50]$
• $y$-axis: rotation parameter $k \in [0, 1]$
• colour: value of $|\delta|$ (dark = zero, bright = maximum)

The light and dark bands are continuous — zero locations migrate smoothly with $k$.
White dotted vertical lines mark the known zeros of $\zeta$ ($k=1$, upper edge of the heatmap).
The transition from the $\beta$-structure (lower edge) to the $\zeta$-structure (upper edge)
is a visual demonstration of the rotation identity:

$$\boxed{\delta(k, s, \infty) = \beta(s)\cdot\cos\!\left(\frac{\pi k}{2}\right) + \frac{\eta(s)}{2^s}\cdot\sin\!\left(\frac{\pi k}{2}\right)}$$

def delta_k(k, s, n_terms=1000):
    """
    Zobecněná δ-funkce s parametrem rotace k ∈ [0,1]:
      k=0   → β(s)             — Dirichletova beta, Catalan G
      k=0.5 → (β+η/2^s)·√2/2  — symetrická kombinace obou
      k=1   → ζ(s)             — Riemannova zeta, π²
    δ(k,s) = Σ_n 2·(n^s-1)/n^s · sin(π(k-n)/2) · (π/2)^(1-s)
    """
    total = complex(0, 0)
    C = (np.pi / 2) ** (1 - s.real)
    for n in range(1, n_terms + 1):
        total += 2 * ((n**s - 1) / n**s) * np.sin(np.pi * (k - n) / 2) * np.pi**(s-1) * C
    return total

# Graf 1: Tři křivky k=0, 0.5, 1
t_range = np.linspace(0, 50, 800)
k_values = [0, 0.5, 1.0]
labels   = ['k=0 → β(s)', 'k=0.5 → mix', 'k=1 → ζ(s)']
colors   = [pastel_coral, pastel_orange, pastel_blue]

results  = {k: [delta_k(k, complex(0.5, t), 800) for t in t_range] for k in k_values}
zeta_std = [riemann_zeta(complex(0.5, t)) for t in t_range]

fig, axes = plt.subplots(3, 1, figsize=(16, 12))
for ax, component, ylabel, title in zip(
        axes,
        [abs, lambda z: z.real, lambda z: z.imag],
        ['|δ|', 'Re δ', 'Im δ'],
        ['Modulo na kritické přímce', 'Re', 'Im']):
    ax.plot(t_range, [component(z) for z in zeta_std],
            color='black', linewidth=1.5, label='ζ', alpha=0.6)
    for k, label, color in zip(k_values, labels, colors):
        ax.plot(t_range, [component(v) for v in results[k]],
                color=color, linewidth=1.8, label=f'δ({label})', linestyle='--', alpha=0.9)
    for zt in zeros_t:
        ax.axvline(x=zt, color='green', linestyle=':', linewidth=0.6, alpha=0.5)
    ax.axhline(y=0, color='black', linewidth=0.5)
    ax.set_ylabel(ylabel, fontsize=12)
    ax.set_title(title, fontsize=13, fontweight='bold')
    ax.legend(fontsize=10, ncol=4); ax.grid(True, alpha=0.3)
    ax.set_xlim(0, 50)
plt.tight_layout()
plt.show()

# Graf 2: Heatmapa |δ(k, 0.5+it)| v prostoru (k, t)
k_range  = np.linspace(0, 1, 60)
t_range2 = np.linspace(0, 50, 400)
heatmap  = np.array([[abs(delta_k(k, complex(0.5, t), 500))
                      for t in t_range2] for k in k_range])

fig, ax = plt.subplots(figsize=(18, 7))
im = ax.imshow(heatmap, aspect='auto', origin='lower',
               extent=[0, 50, 0, 1], cmap='plasma', vmin=0, vmax=3)
plt.colorbar(im, ax=ax, label='|δ(k, 0.5+it)|')
for zt in zeros_t:
    ax.axvline(x=zt, color='white', linestyle=':', linewidth=0.8, alpha=0.7)
ax.axhline(y=0,   color='cyan',   linewidth=1.5, label='k=0 → β(s)')
ax.axhline(y=0.5, color='yellow', linewidth=1.5, label='k=0.5 → mix')
ax.axhline(y=1.0, color='lime',   linewidth=1.5, label='k=1 → ζ(s)')
ax.set_xlabel('t', fontsize=13); ax.set_ylabel('k (parametr rotace)', fontsize=13)
ax.set_title('Heatmapa |δ(k, 0.5+it)|: rotace od β (k=0) k ζ (k=1)',
             fontsize=14, fontweight='bold')
ax.set_yticks([0, 0.25, 0.5, 0.75, 1.0])
ax.set_yticklabels(['0 (β)', '0.25', '0.5 (mix)', '0.75', '1 (ζ)'])
ax.legend(fontsize=11, loc='upper right')
plt.tight_layout()
plt.show()

24. Why Even Numbers Have Modulus $> 0$ Off the Critical Line

Minimum Exactly at $\sigma = \frac{1}{2}$

The unnormalised $\delta_{\text{even}}$ (without division by $(1-2^{1-s})$) exhibits a
global minimum of the modulus on the critical line $\sigma = \frac{1}{2}$,
even without any normalisation — arising purely from the structure of the even sum.

At the first zero of $\zeta$ ($t \approx 14.135$):

| $\sigma$ | $|\delta_\text{even}|$, $N=3000$ | |----------|----------------------------------| | 0.1 | 0.298 | | 0.3 | 0.176 | | 0.5 | 0.007 ← minimum | | 0.7 | 0.230 | | 0.9 | 0.493 |

The value decreases monotonically towards $\sigma = 0.5$ and grows symmetrically on both sides.
As $N$ increases, the minimum at $\sigma = \frac{1}{2}$ deepens — converging to zero.

Why: Symmetry of the Weights $|n^{-s}|$ on the Critical Line

Each term in the sum carries the weight $\left(1 - n^{-s}\right)$, where $|n^{-s}| = n^{-\sigma}$:

$$\left|n^{-(\sigma+it)}\right| = n^{-\sigma} \cdot \underbrace{|e^{-it\ln n}|}_{=1}= n^{-\sigma}$$

On the critical line $\sigma = \frac{1}{2}$ we have $|n^{-s}| = n^{-1/2} = \frac{1}{\sqrt{n}}$ — the weights
are the geometric mean of the values at $\sigma = 0$ and $\sigma = 1$.

This is precisely the value around which the functional equation is symmetric:

$$\zeta(s) = 2^s \pi^{s-1} \sin\!\left(\frac{\pi s}{2}\right)\Gamma(1-s)\,\zeta(1-s)$$

The symmetry $s \leftrightarrow 1-\bar{s}$ is centred on $\sigma = \frac{1}{2}$.
For $\sigma \neq \frac{1}{2}$ this symmetry of the weights $n^{-\sigma}$ and $n^{-(1-\sigma)}$ is broken
— the oscillatory cancellation of the phase contributions $e^{-it\ln n}$ is imperfect and $|\delta_\text{even}| > 0$.

Algebraic Cause: The Zeros of $\delta_\text{even}$ Are Not the Zeros of $\zeta$

In the limit $N \to \infty$ (with regularisation):

$$\delta_\text{even}(s) = \pi^{s-1}\!\left[\eta(s) - 2^{s-1}\right]$$

$\delta_\text{even}(s) = 0$ occurs when $\eta(s) = 2^{s-1}$ — not when $\eta(s) = 0$.

At the zeros of $\zeta$ we have $\eta(s_0) = 0$, but $2^{s_0-1} \neq 0$ (on the critical line $|2^{s-1}| = \frac{1}{\sqrt{2}}$):

$$\boxed{\delta_\text{even}(s_0) = \pi^{s_0-1}\!\left[0 - \frac{1}{\sqrt{2}}\right] \neq 0}$$

Only the normalisation $\frac{1}{1-2^{1-s}}$ removes this additive constant and restores
agreement with the zeros of $\zeta$. Without normalisation, the zeros of $\delta_\text{even}$ lie at different locations from those of $\zeta$
— and off the critical line, $|\delta_\text{even}|$ grows larger the further $\sigma$ is from $\frac{1}{2}$.

#HERE IS A ADDITIONAL FUNCTION FOR COMPUTING ZEROS
import numpy as np
from scipy.optimize import minimize_scalar
from scipy.signal import find_peaks

def modified_zeta_v3(s, n_terms=2000):
    """
    δ-funkce: přesná reprezentace ζ(s) přes sudá čísla
    δ(s) = Σ_{n sudé} 2·(n^s-1)/n^s · sin(π(0-n)/2) · (π/2)^(1-s) · 1/(1-2^(1-s))
    Algebraicky: δ(s) = ζ(s)
    """
    total = complex(0, 0)
    C = (np.pi / 2) ** (1 - s.real)     # přesný koeficient C(s) = (π/2)^(1-s)

    for n in range(2, n_terms + 1):
        if n % 2 == 0:
            total += ((n**s/8**s - 1) / n**s) * np.sin(np.pi * (0 - n) / 8) * np.pi**(s-1) * C

    total = total

    return total

def find_zeros(max_t=1000, resolution=0.02, threshold=0.01):
    """Find zeros on critical line s = 0.8 + it"""
    
    # Scan
    t_scan = np.arange(0, max_t, resolution)
    abs_values = []
    
    print(f"Scanning t=0 to {max_t}...")
    for i, t in enumerate(t_scan):
        if i % 5000 == 0:
            print(f"  Progress: {100*i/len(t_scan):.1f}%")
        s = complex(0.8, t)
        z = modified_zeta_v3(s, n_terms=2000)
        abs_values.append(abs(z))
    
    abs_values = np.array(abs_values)
    
    # Find minima
    peaks_inv, _ = find_peaks(-abs_values, prominence=0.001, distance=int(0.3/resolution))
    
    # Refine
    zeros = []
    for peak_idx in peaks_inv:
        t_approx = t_scan[peak_idx]
        
        if any(abs(t_approx - z) < 0.2 for z in zeros):
            continue
        
        def objective(t):
            return abs(modified_zeta_v3(complex(0.8, t), n_terms=8000))
        
        try:
            result = minimize_scalar(objective, 
                                    bounds=(t_approx - 0.8, t_approx + 0.8),
                                    method='bounded')
            if result.fun < threshold:
                zeros.append(result.x)
        except:
            pass
    
    return sorted(zeros)

# Find zeros
zeros = find_zeros(max_t=2300, resolution=0.02, threshold=0.005)


print(f"Found {len(zeros)} zeros")

print("modified_zeta_v3_zeros = [")
for i in range(0, len(zeros), 5):
    chunk = zeros[i:i+5]
    if i + 5 < len(zeros):
        print("    " + ", ".join(f"{z:.10f}" for z in chunk) + ",")
    else:
        print("    " + ", ".join(f"{z:.10f}" for z in chunk))
print("]")