```#!/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.
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
```