#!/usr/local/bin/python3

from math import exp, cos
from numpy import pi

def f(x):
    return exp(-x * x)

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

# This class is used internally to hold stack data for adaptive_simpson_integration.
class Adaptive_vector:
    def __init__(self, a, h, fa, faph, fap2h, area):
        self.a, self.h, self.fa, self.faph, self.fap2h, self.area = a, h, fa, faph, fap2h, area
    
def adaptive_simpson_integration(f, lower_limit, upper_limit, epsilon = 1.0e-7, stack_max = 100):
    delta, integral, middle, vs = upper_limit - lower_limit, 0.0, (lower_limit + upper_limit) / 2.0, []
    h = delta / 2.0
    f_lower_limit, f_upper_limit, f_middle = f(lower_limit), f(upper_limit), f(middle)
    area = (f_lower_limit + 4 * f_middle + f_upper_limit) * h / 3.0
    vs.append(Adaptive_vector(lower_limit, h, f_lower_limit, f_middle, f_upper_limit, area))
    while True:
        # Save repeated indexing
        top_v = vs.pop()
        h = top_v.h / 2.0
        # Two new points needed, f(middle of left half) and f(middle of right half)
        f_left_middle = f(top_v.a + h)
        area_left = (top_v.fa + 4 * f_left_middle + top_v.faph) * h / 3.0
        f_right_middle = f(top_v.a + 3.0 * h)
        area_right = (top_v.faph + 4 * f_right_middle + top_v.fap2h) * h / 3.0
        if abs(area_left + area_right - top_v.area) < 60 * epsilon * h / delta:
            # Function is smooth enough to accept left and right half areas.
            integral += area_left + area_right + (area_left + area_right - top_v.area) / 15.0
            if not vs:
                break
        else:
            # Function is NOT smooth enough to accept left and right half areas.
            # Make 2 new "half" vectors and update the stack.
            a_v1 = Adaptive_vector(top_v.a, h, top_v.fa, f_left_middle, top_v.faph, area_left)
            a_v2 = Adaptive_vector(top_v.a + 2.0 * h, h, top_v.faph, f_right_middle, top_v.fap2h, area_right)
            vs.append(a_v1)
            vs.append(a_v2)
            if len(vs) >= stack_max:
                raise RuntimeError('Stack exceeded', stack_max)

    return integral

if __name__ == "__main__":
    import sys
    try:
        integral = 2.0 * adaptive_simpson_integration(f, 0.0, 8.0, epsilon = 1.0e-13)
        print('Integral**2 = {0:15.13f}, Exact = {1:15.13f}'.format(integral * integral, pi))
        integral = adaptive_simpson_integration(gamma, 0.0, 90.0, epsilon = 1.0e-13)
        print('\n5! = Gamma(6) =~ {0:14.12f}'.format(integral))
        integral = adaptive_simpson_integration(cos, 0.0, pi / 2.0, epsilon = 1.0e-13)
        print('\nCosine from 0 to pi / 2 = {0:14.12f}'.format(integral))
    except RuntimeError as e:
        print(e)
# End of Code

# Test Ouput:
# Integral**2 = 3.1415926535898, Exact = 3.1415926535898

# 5! = Gamma(6) =~ 120.000000000000

# Cosine from 0 to pi / 2 = 1.000000000000