#
# PYTHON • MyCompiler
# eulermaclaurin
# https://[Log in to view URL]
# Euler-Maclaurin
# 
# 2026 Maio Dom, 31 — Junho Ter, 2
# by 
# j
# 

import matplotlib.pyplot as plt

import numpy as np

from numpy import pi
from numpy import angle

from scipy.special import zeta
from scipy.special import bernoulli

import math
import colorsys

#
# W A R N I N G 
#
# Para obter gráficos de melhor qualidade
# é necessário configurar a variável 
# MAXRESOLUTION para valores tal qual 
# 512 ou 1024; no entanto haverá a 
# necessidade de mudar a ambiência de 
# execução deste código python para 
# evitar o estouro de tempo.
#
# To obtain higher-quality graphics,
# you need to set the MAXRESOLUTION 
# variable to values such as 512 or
# 1024; however, you will need to 
# change the execution environment 
# of this Python code to avoid a 
# timeout.
#

MAXRESOLUTION = 72

def eulermaclaurin(prt,s):

    #
    # Definições
    # 

    N = 40
    M = 36

    ber = bernoulli(2*M)
    
    #
    # Cálculo Primeiro Termo
    #

    s1 = 0
    for n in range(1,N):
        s1 += n**(-s)

    #
    # Cálculo Segundo Termo
    #

    s2 = (N**(1-s))/(s-1)

    #
    # Cálculo Terceiro Termo
    #

    s3 = (1/2)*N**(-s)

    #
    # Cálculo Quarto Termo
    #

    s4 = 0
    for k in range(1,M+1):
        f = ber[2*k]/math.factorial(2*k)
        num = 1
        max = 2*k - 2 + 1
        for j in range(0,max):
            num = num*(s + j)
        den = N**(s+2*k-1)
        rz = num/den
        s4 += f*rz
        
    #
    # Cálculo Geral | Função zeta
    # 
    
    zz = s1 + s2 + s3 + s4

    zr = zeta(s)

    #
    # Resultados
    # 

    if (prt==1):
        
        print('s1   = '+f'{s1:.6f}')
        print('s2   = '+f'{s2:.6f}')
        print('s3   = '+f'{s3:.6f}')
        print('s4   = '+f'{s4:.6f}')
        print('zeta = '+f'{zz:.12f}')
        print('zeta = '+f'{zr:.12f}')

    #
    # Retorno
    # 

    return zz

def Graphriemann():
    
    #
    # Parametrizações
    #
    
    plt.style.use('dark_background')
    
    xRESOLUTION = MAXRESOLUTION
    yRESOLUTION = MAXRESOLUTION
    
    H   = np.zeros((xRESOLUTION, yRESOLUTION))
    S   = np.zeros((xRESOLUTION, yRESOLUTION))
    L   = np.zeros((xRESOLUTION, yRESOLUTION))
    RGB = np.zeros((xRESOLUTION, yRESOLUTION, 3))
    
    xMIN = -5.00
    yMIN = -5.00
    xMAX = +5.00
    yMAX = +5.00

    #
    # Preparação
    #
    
    deltaX = (xMAX - xMIN) / xRESOLUTION
    deltaY = (yMAX - yMIN) / yRESOLUTION

    #
    # Varredura Principal
    #
    
    for idx in range(xRESOLUTION):
        
        x = xMIN + idx * deltaX
        
        for idy in range(yRESOLUTION):
            
            y = yMIN + idy * deltaY
            
            #
            # Cálculo da Função Zeta de Riemann
            #
            
            # s = complex(zeta(complex(x, y)))
            
            try:
                
                # mp.zeta suporta números complexos em todo o plano
                #z = complex(mp.zeta(s))
                
                s = complex(eulermaclaurin(0,complex(x, y)))
                # s = complex(zeta(complex(x, y)))
                
            except:
                
                s = 0j
                
            #
            # Cálculo parâmetros cores HLS
            #
            
            H[idx, idy] = (angle(s) + pi) / (2 * pi)
            L[idx, idy] = 1.0 - 1.0 / (1.0 + abs(s)**0.5)
            S[idx, idy] = 1.0

    #
    # Converte cada pixel individualmente 
    # de HLS para RGB 
    # usando a biblioteca colorsys
    #
    
    for i in range(xRESOLUTION):
        
        for j in range(yRESOLUTION):
        
            RGB[i, j] = colorsys.hls_to_rgb(H[i, j], L[i, j], S[i, j]) 

    #
    # Gera o gráfico
    # 
    
    plt.imshow(RGB.swapaxes(0, 1), extent=[xMIN, xMAX, yMIN, yMAX], origin='lower')

    plt.title("Riemann Zeta function $ζ(s)$")
    plt.xlabel("Re(s)")
    plt.ylabel("Im(s)")
    
    plt.axis('on')
    plt.show()

eulermaclaurin(1,complex(1/2,135))
Graphriemann()