#!/usr/local/bin/python3

from math import exp, pi, sqrt

import sys
sys.path.append('/Users/jim/Python_3/numerical_libraries')

from diff_and_int import adaptive_simpson_integration, Adaptive_vector

def gamma(t, z = 1.75):
    return t**(z - 1.0) * exp(-t)

# From Artin's book
x = 4.75
print('x =', x)

x_sq = x * x

# Eq. 5.20
mu = ((((1.0 / (1260.0 * x_sq)) - 1.0 / 360.0) / x_sq) + 1.0 / 12.0) / x
print('mu = {0:.6f}'.format(mu))

# Stirling's Formula, Eq. 3.19 (first equation)
g_x = sqrt(2 * pi) * x**(x - 1.0 / 2.0) * exp(mu - x)
print('gamma({0:.2f}) = {1:.6f}'.format(x, g_x))

arg = 1.75
result = g_x / 3.75 / 2.75 / arg
print('\nFrom Artin   : gamma({0:.2f}) = {1:.7f}'.format(arg, result))

integral = adaptive_simpson_integration(f = gamma, lower_limit = 0.0, upper_limit = 20.0, epsilon = 1.0e-6)
print('From integral: gamma({0:.2f}) = {1:.7f}'.format(arg, integral))

# Output:

# x = 4.75
# mu = 0.017518
# gamma(4.75) = 16.586207

# From Artin   : gamma(1.75) = 0.9190625
# From integral: gamma(1.75) = 0.9190625