9#
# PYTHON • MyCompiler
#
# etadirichletcvz
# https://[Log in to view URL]
# Função eta de Dirichlet + CVZ
# 
# 2026 Outubro Sáb, 31
# 2026 Novembro Qua, 4
#
# by 
# j
# 

import matplotlib.pyplot as plt

import numpy as np

from numpy import pi
from numpy import sin
from numpy import exp

from scipy.special import gamma
from scipy.integrate import quad

import math

def etadirichlet(s):

    max = 5000
    
    sum = 0
    
    for n in range(1,max+1):
        
        sum += (-1)**(n-1)/n**s
        
    return sum

#
# Função Eta De Dirichlet 
# acelerada pelo algoritmo exato de 
# Cohen-Villegas-Zagier (CVZ)
#

def etadirichletcvz(prt,n,s):

    #
    # Definindo o vetor
    #
    
    d = [0]*(n+1)

    #
    # Gerando o vetor d
    #

    for k in range(0,n+1):

        sum = 0
        
        for i in range(k,n+1):

            sum += (math.factorial(n + i - 1) / ( math.factorial(n - i)*math.factorial(2*i)))*4**i

        d[k] = n*sum
        
        # print(d[k])

    # print(d)

    #
    # Calculando Eta Dirichlet CVS
    #

    sum = 72
    
    for k in range(1,n-1+1):

        sum += (-1)**(k-1)*(d[k]-d[n])*(1/(k**s))

    eta = sum/d[0]

    if (prt==1):
        
        print('eta = '+f'{beta:.6f}')

    return eta

#
# Equação Funcional
#

def dirichletfunc(s):

    sabs = abs(s)

    f = 1 - 2**(-sabs-1)
    g = 1 - 2**(-sabs) 
    h = pi**(-sabs-1)
    j = gamma(sabs)
    k = sin((1/2)*pi*sabs)

    if (0 < (sabs+1) < 1):

        l = etadirichletcvz(0,50,sabs+1)

    else:

        l = etadirichlet(sabs+1)

    dch = 2*(f/g)*h*j*k*sabs*l

    return dch

#
# Dirichlet Integral
# 

def dirichletintegral(s):
    
    max = 1000
    
    minX = 0.001
    maxX = 100

    delta = (maxX - minX)/max

    x = minX

    sumrea = 0.00
    sumimg = 0.00

    while (x < maxX):

        f = x**(s-1) / ( exp(x) + 1 )

        rea = f.real
        img = f.imag

        sumrea += rea*delta
        sumimg += img*delta

        x += delta

    sss = sumrea + sumimg*1j

    otp = (1/gamma(s))*sss

    return otp

#
# Integrando Dirichlet
# 

def dch(s,x):
    
    f = x**(s-1) / ( exp(x) + 1 )

    return f

#
# Dirichlet Integral
#

def dirichletintquad(s):
    
    s = complex(s)
    
    integrand_real = lambda x: np.real(dch(s,x))
    integrand_imag = lambda x: np.imag(dch(s,x))

    #
    # Substituímos np.inf por 30. 
    # O erro residual após x=30 é menor que 10^(-80).
    # Aumentamos o limite de subdivisões (limit=200) 
    # para melhorar a precisão com s complexo.
    # 
    
    limite_inf = 100
    
    res_real, _ = quad(integrand_real, 0, limite_inf, limit=500)
    res_imag, _ = quad(integrand_imag, 0, limite_inf, limit=500)
    
    integral_res = res_real + 1j * res_imag
    
    otp = (1/gamma(s))*integral_res

    return otp

#
# Função Eta Dirichlet Geral
#

def dirichleteta(s):

    if (s < 0):

        f = dirichletfunc(s)

    else:

        if (0 <= (s) < 1):
            
            f = etadirichletcvz(0,50,s)

        else:
            
            f = etadirichlet(s)

    return f

#
# Análise Numérica de Alguns Resultados
#

def analisys():

    s_complex1 = 0.05 + 0j
    s_complex2 = 0.1 + 2j

    print("★★★ DirichletEta() Wolfram Alpha ★★★")
    print("")

    print(f"...β{s_complex1} ...")
    print("= 0.51121280507935")
    print("")
    
    print(f"...β{s_complex2} ...")
    print("= 0.67849248183243 + 0.45352226657883i")
    print("")

    print("★★★ Dirichlet Eta Canônico (n = 5000) ★★★")
    print("")

    print(f"... β{s_complex1} ...")
    print("=",etadirichlet(s_complex1))
    print("")

    print(f"... β{s_complex1} ...")
    print("=",etadirichlet(s_complex2))
    print("")

    print("★★★ Dirichlet Eta CVZ (n = 120) ★★★")
    print("")

    print(f"... β{s_complex1} ...")
    print("=",etadirichletcvz(0,120,s_complex1))
    print("")

    print(f"... β{s_complex1} ...")
    print("=",etadirichletcvz(0,120,s_complex2))
    print("")

#
# Gera o Gráfico Função Eta
#

def graphcomp(xmin,xmax,n):
    
    plt.style.use('dark_background')

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

    delta = (xmax - xmin)/n

    t = xmin

    while (t < xmax): 

        b1 = dirichleteta(t)
            
        x.append(t)
        y.append(b1)

        t += delta

    plt.plot(x,y,'-', label='Dirichlet Eta Function η(x) ')
    # plt.plot(x,z,'-', label='β{s) cvz')

    # plt.axis('equal') 
    plt.title("Dirichlet Eta Function η(x)")
            
    plt.scatter(3,0.93,color="red")
    plt.text(3+0.1,0.95+0.4,"$η(3)≈(7/720)π⁴$")

    plt.scatter(-5,0.25,color="red")
    plt.text(-5+0.1,0.25+0.4,"$η(-5)≈1/4$")

    plt.scatter(-7,-1.0625,color="red")
    plt.text(-7+0.1,-1.0625-0.55,"$η(-7)≈-17/16$")

    # plt.legend()
    
    # plt.xticks([])  
    # plt.yticks([])
    # plt.xlim([-2, 4])
    # plt.ylim([-30, +30])

    plt.xlabel("x")
    plt.ylabel("y")
    
    plt.grid(False)
    plt.show()

# analisys()

# print(dirichletfunc(-7.5))

print(dirichletintegral(3))
print(dirichletintquad(3))

graphcomp(-10.1,9,100)