#
# PYTHON • MyCompiler
#
# betadirichletcvz
# https://[Log in to view URL]
# Função Beta de Dirichlet + CVZ
# 
# 2026 Outubro Ter, 20
# 2026 Outubro Sex, 30
#
# by 
# j
# 

import matplotlib.pyplot as plt

import numpy as np

from numpy import pi
from numpy import angle
from numpy import cos

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

import cmath
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 
# 256 ou 512; 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 256 or
# 512; however, you will need to 
# change the execution environment 
# of this Python code to avoid a 
# timeout.
#
# P U R P O S E
# 
# O propósito deste código é apresentar
# o gráfico complexo da Função Beta de
# Dirichlet quando usamos a Função 
# Canônica para o cálculo em Re(s) > 0
# que gera erros nos primeiros valores
# apresentado um gráfico com uma 
# incômoda faixa & também, quando usamos
# a Função de aceleração 
# Cohen-Villegas-Zagier (CVZ) que conduz
# a um gráfico perfeito. Os gráficos 
# podem ser gerados mediante os 
# comandos ...
#
# Função Canônica
# menu(3,0)
#
# Função com CVZ
# menu(3,1)
# 
# The purpose of this code is to display
# the complex graph of the Dirichlet Beta 
# Function when using the Canonical 
# Function for the calculation in 
# Re(s) > 0, which generates errors
# in the first few values, resulting 
# in a graph with an annoying band; and 
# also, when using the (CVZ), which 
# leads to a perfect graph. The graphs 
# can be generated using the following 
# commands...
# 
# Canonical Function
# menu(3,0)
# 
# Function with CVZ
# menu(3,1)
# 

MAXRESOLUTION = 60

#
# Função Beta De Dirichlet Canônica
# Re(s) > 0
#

def betadirichletR(prt,s):

    max = 1000

    sum = 0

    for k in range(0, max + 1):
        
       sum += (-1)**k / (2*k + 1)**s

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

    return sum

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

def betadirichletRcvz(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 Beta Dirichlet CVS
    #

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

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

    beta = sum/d[0]

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

    return beta

#
# Extensão Analítica da Função
# Beta de Dirichlet
# Re(s) < 0
#

def betadirichletL(prt,s):
    
    #
    # Cálculo Termos
    #

    p1 = ( (1/2)*pi )**(s-1)
    p2 = gamma(1-s)
    p3 = cos((1/2)*pi*s)
    p4 = betadirichletR(0,1-s)

    #
    # Cálculo Beta
    #

    bb = p1*p2*p3*p4
    
    #
    # Resultados
    # 

    if (prt==1):
        
        print('p1   = '+f'{p1:.6f}')
        print('s2   = '+f'{p2:.6f}')
        print('s3   = '+f'{p3:.6f}')
        print('p4   = '+f'{p4:.6f}')
        print('beta = '+f'{bb:.12f}')

    #
    # Retorno
    # 

    return bb

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

def analisys():

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

    print("★★★ DirichletBeta() Wolfram Alpha ★★★")
    print("")

    print(f"...β{s_complex1} ...")
    print("= 0.51928929459480")
    print("")
    
    print(f"...β{s_complex2} ...")
    print("= 1.0660051190064 + 0.5900188380572")
    print("")

    print("★★★ Dirichlet Beta Canônico ★★★")
    print("")

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

    print(f"... β{s_complex2} ...")
    print("=",betadirichletR(0,s_complex2))
    print("")

    print("★★★ Dirichlet Beta + cvz (n=116) ★★★")
    print("")

    print(f"... β{s_complex1} ...")
    print("=",betadirichletRcvz(0,116,s_complex1))
    print("")
    
    print(f"... β{s_complex2} ...")
    print("=",betadirichletRcvz(0,116,s_complex2))
    print("")

#
# Verifica Ponto Mínimo
#

def checkpoint():

    eps = 0.001
    min = 1

    n = 10000
    xMIN = 0.825
    xMAX = 0.827
    
    deltaX = (xMAX - xMIN)/n

    f = True

    x = xMIN
    while ( (x < xMAX) and f ):
        
        b1 = betadirichletR(0,x)
        b2 = betadirichletRcvz(0,25,x)
        
        d = abs(b1 - b2)
        
        if d < min:
        	
        	min = d
        	xk = x

        if d < eps:

            f = False

        x += deltaX

    if not(f):

        print("★★★ Conclusão ★★★")
        print("")
        print(f"A partir do ponto x={x}")
        print("o erro absoluto entre β(x) canônico")
        print(f"e β(x) cvz é menor que {eps}")
        print("")
        print(f"β({x:.5f}) Canônico = {betadirichletR(0,x):.3f}")
        print(f"β({x:.5f}) cvz      = {betadirichletRcvz(0,25,x):.3f}")
        print("")
              
    else:

        print("Não Encontrado!")
    
    # print("xk =",xk)
    # print("min =",min)
    
#
# Gera o Gráfico Comparativo
#

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

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

    delta = (xmax - xmin)/n

    t = xmin

    while (t < xmax): 

        b1 = betadirichletR(0,t)
        b2 = betadirichletRcvz(0,25,t)
            
        x.append(t)
        y.append(b1)
        z.append(b2)

        t += delta

    plt.plot(x,y,'-', label='β(s) canônico')
    plt.plot(x,z,'-', label='β{s) cvz')

    # plt.axis('equal') 
    plt.title("β(x) canônico × β(x) cvz")
            
    plt.scatter(0.05,0.86,color="red")
    plt.text(0.05+0.05,0.86+0.00,"$β(0.05)≈0.86$")

    plt.scatter(0.05,0.52,color="red")
    plt.text(0.05+0.05,0.52+0.00,"$β(0.05)=0.52$")
    
    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()
        
def graph(cvz):
    
    #
    # 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 = -4.00
    yMIN = -4.00
    xMAX = +4.00
    yMAX = +4.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 Beta de Dirichlet
            #
            
            try:

                if x <= 0:
                    
                   s = complex(betadirichletL(0,complex(x, y)))

                else:

                   match cvz:

                       case 0:
                     
                         s = complex(betadirichletR(0,complex(x,y)))

                         smode = ""

                       case 1:

                         s = complex(betadirichletRcvz(0,25,complex(x,y)))

                         smode =" + cvz"
                
            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("Dirichlet Beta Function $β(s)$"+smode)
    plt.xlabel("Re(s)")
    plt.ylabel("Im(s)")
    
    plt.axis('on')
    plt.show()

#
# Menu Geral
#

def menu(opt,g):

    match opt:

        #
        # Análises Comparativas
        # 

        case 0:
            
            analisys()

        #
        # Ponto Crítico
        #

        case 1:

            checkpoint()

        #
        # Gráfico Comparativo
        #

        case 2:

            graphcomp(0,0.001,1,100)

        #
        # Gráficos Complexos
        # 
        
        case 3:

            if (g == 0):
                
                graph(0)

            if (g == 1):
                
                graph(1)

menu(0,1)