#
# PYTHON • MyCompiler
# gammainverse
# https://[Log in to view URL]
# Função Gamma Inversa
# 
# 2026 Maio Seg, 25 — Qua, 27
# by 
# j
# 

import matplotlib.pyplot as plt
import numpy as np
from numpy import sqrt
from numpy import exp
from numpy import log
from numpy import pi
from scipy.special import gamma
from scipy.special import lambertw

import math

def pyg(m,z):

    max = 3000

    psi = 0.00
    
    for k in range(0,max+1):
        
        if ( (z+k) != 0 ):
            
            psi += 1/((z + k)**(m+1))
                
    psi = ((-1)**(m+1))*math.factorial(m)*psi

    return psi

def lacomo_boyd(b,x):
  max = 1000
  e = exp(1)
  d = sqrt(2)
  P = sqrt(1+e*x)
  if ((b==0) and (x>-(1/e)) and (x<0)):
    w =(e*x*log(1+P))/(1+e*x+P)
  if ((b==0) and (x>0) and (x<e)):
    w = x/e
  if ((b==0) and (x>e)):
    w = log(x) - log(log(x))
  if ((b==-1) and (x>-(1/e)) and (x<-1/4)):
    w = -1 - d*P
  if ((b==-1) and (x>-(1/4)) and (x<0)):
    w = log(-x)-log(-log(-x))
  for i in range(0,max+1):
    w = (w/(1+w))*(1+log(x/w))
  return w

def W(x):
    xmin = - 1/exp(1)
    if (x > xmin):
        r = lacomo_boyd(0,x)
    else:
        r = lambertw(x)
    return r

def gammainverse(x):
    e = 1/exp(1)
    b = x/(sqrt(2*pi))
    f = log(b)/W(e*log(b)) + 1/2
    return f.real

def newtonraphson(b,x):

    zmin = 1.4616321449683623
    
    max = 500
    
    eps = 1e-12
    
    f = True
    
    if (b == 0):
        
        zn = gammainverse(x)
        sinal = +1
      
    else:
     
        zn = 1.0 / x if x > 2 else 0.5
        sinal = -1
     
    n = 0
    
    while ( ( n < max ) and f ):
        
        zp = zn + sinal*((1-x/gamma(zn))/pyg(0,zn))
        
        if (b == -1):
            
            if ( ( zp <= 0 ) or (zp >= zmin) ):
                
                zp = zn / 2.0 
       
        if ( abs(x-gamma(zp)) < eps ):
            
            f = False
            
        else:
            
            zn = zp
            n += 1
            
    return zp
            
def graf():

    plt.style.use('dark_background')

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

    xmin = 1
    xmax = 1000
    n = 100

    delta = (xmax - xmin)/n

    t = xmin

    while (t < xmax): 
        
        x.append(t)
        y.append(gammainverse(t))

        t += delta

    plt.plot(x,y,'-')
        
    tit = 'Inverse Gamma Function'

    # 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()

def graph():
    
    plt.style.use('dark_background')

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

    xmin = 0.009
    xmax = 6
    n = 50

    delta = (xmax - xmin)/n

    t = xmin

    while (t < xmax): 
        
        x.append(gamma(t))
        y.append(t)

        t += delta

    plt.plot(x,y,'-')
        
    tit = 'Inverse Gamma Function'

    # plt.axis('equal') 
    plt.title(tit)
    
    plt.scatter(24,5,color="red")
    plt.text(24+0.35,5-0.45,"$Γ^{-1}(24)=5$")

    plt.scatter(24,0.1,color="red")
    plt.text(24+0.35,0.04075+0.35,"$Γ^{-1}(24)≈0.0407532$")

    plt.scatter(0.8856,1.4616,color="yellow")
    plt.text(0.8856+3,1.4616+0.0,"Junction point of the branches")

    plt.text(0.8856+3,1.4616-0.35,"(0.8856031944108889,1.4616321449683623)")
        
    # plt.legend()
    # plt.xticks([])  
    # plt.yticks([])
    # plt.xlim([-2, 4])
    # plt.ylim([-pico, pico])
    plt.grid(False)
    plt.show()

def view(x):

    print("... [ Ramo Principal ] Via Stirling ...")
    print("")
    print("GammaInverse(",x,") =",gammainverse(x))
    
    print("")
    print("... [ Ramo Primciapal ] Via Newton-Raphson ...")
    print("")
    print("GammaInverse(",x,") =",newtonraphson(0,x))

    print("")
    print("... [ Ramo Secundário ] Via Newton-Raphson ...")
    print("")
    print("GammaInverse(",x,") =",newtonraphson(-1,x))

    print("")
    print('*** =',gamma(newtonraphson(-1,x)))
    
    graph()

view(24)

#gi = 0.8856031944108889
#gi = 0.90
#print('gi == ',gi)
#g=newtonraphson(0,gi)
#print('Inverse = ',g)
#print('check = ',gamma(g))

#print('Referencia = 1.4616321449683623')