#
# PYTHON • MyCompiler
# riemannsiegel
# https://[Log in to view URL]
# Riemann-Siegel
# 
# 2026 Maio Qui, 28 — Sáb, 30
# by 
# j
# 

import matplotlib.pyplot as plt
import numpy as np

from numpy import sqrt
from numpy import cos
from numpy import exp
from numpy import log
from numpy import pi
from numpy import floor

from scipy.special import gamma
from scipy.special import loggamma
from scipy.special import zeta

import math

def c0(p):
    
    den = cos(2 * pi * p)
    
    if (abs(den) < 1e-11):

        c = 0

    else:
        
        num = cos(2 * pi * (p**2 - p - 1/16))
        
        c = num / den

    return c

def c1(p):

    num = 136 - 495*p + 405*p**2 - 90*p**3

    den = 1080

    c = num / den

    return c

def c2(p):
    
    num = 135 - 3960*p + 16200*p**2 - 21600*p**3 + 8640*p**4
    
    den = 466560
    
    c = num / den
    
    return c

def c3(p):
    
    x = 2*p - 1

    A = 0.04342
    B = -0.3211
    D = 0.7412

    c = (1/(720*8))*x*(x**2-1)*(A*x**4+B*x**2+D)

    return c

def c4(p):

    x = 2*p - 1

    A = 0.00121
    B = -0.145
    D = 0.521
    E = -0.412

    c = (1/(40320*16))*(x**2-1)*(A*x**6+B*x**4+D*x**2+E)
    
    return c
    
def riemann(prt,t):
    
    #
    # Cálculo do número N de iterações
    # 

    N = int(floor(sqrt(t/(2*pi))))

    #
    # Cálculo do teta
    #

    # teta = (t/2)*log(t/(2*pi))-(t/2)-pi/8 + 1/(48*t) + 7/(5760*t**3)

    teta = loggamma(1/4 + (t/2)*1j).imag - (t/2)*log(pi)
    
    #
    # Estimativa Assintótica Z(t)
    #

    s = 0
    for k in range(1,N+1):
        s += (1/sqrt(k))*cos(teta - t*log(k))

    s = 2*s

    #
    # Cálculo de p
    # 

    p = sqrt(t / (2 * pi)) - N

    #
    # Cálculo de R
    # 

    tau = (2*pi)/t

    f = ((-1)**(N - 1)) * tau**0.25 
    
    # R = f * (1*c0(p) + tau**(1/2)*c1(p) + tau**(1)*c2(p) + tau**(3/2)*c3(p)+tau**(2)*c4(p))

    R = f * (1*c0(p) + tau**(1/2)*c1(p) + tau**(1)*c2(p) )
    
    #
    # Ajuste de Z(t)
    #

    z = s + R
    
    #
    # Cálculo Função zeta
    #
    # z(½ + it) = e^(-iteta)*Z(t)
    #

    zet = exp(-1j*teta)*z

    #
    # Valor Real
    #

    zett = zeta(1/2 + 135*1j)

    #
    # Apresenta Resultado
    #

    if (prt == 1):
        
        print('t         =',t)
        print('N         =',N)
        print('teta      = '+f'{teta:.5f}')
        print('soma      = '+f'{s:.5f}')
        print('p         = '+f'{p:.5f}')
        print('R         = '+f'{R:.5f}')
        print('z         = '+f'{z:.5f}')
        print('zeta      = '+f'{zet:.2f}')
        print('zeta real = '+f'{zett:.5f}')

    return zet
            
def graph():

    plt.style.use('dark_background')

    plt.figure(figsize=(6,6))
    
    x = []
    y = []

    xmin = 0.1
    xmax = 135
    n = 3000

    delta = (xmax - xmin)/n

    t = xmin

    while (t < xmax): 
        
        x.append(riemann(0,t).real)
        y.append(riemann(0,t).imag)

        t += delta

    plt.plot(x,y,'-',color="red")
        
    tit = 'Riemann Zeta Function | s = ½ + it'

    # plt.axis('equal') 
    plt.title(tit)
    # plt.legend()
    # plt.xticks([])  
    # plt.yticks([])
    # plt.xlim([-2, 4])
    # plt.ylim([-pico, pico])
    plt.grid(False)
    plt.show()

riemann(1,135)
graph()