simpson.py

#!/usr/bin/env python3

from math import sqrt, pi

def integrate(f, x, error_margin=0.0001, start_num_seg=10):
    def p(num_seg):
        width = x / num_seg

        p = f(0)
        p += sum(4*f(i*width) for i in range(1, num_seg, 2))
        p += sum(2*f(i*width) for i in range(2, num_seg-1, 2))
        p += f(x)
        p *= width / 3

        return p

    old, new = error_margin + 1, error_margin

    num_seg = start_num_seg
    while abs(new - old) > error_margin:
        old, new = p(num_seg), p(num_seg*2)
        num_seg *= 2

    return new

def t(x, dof):
    def gamma(x):
        if x == 1: return 1
        if x == 1/2: return sqrt(pi)
        return (x-1) * gamma(x-1)
    ret = gamma((dof + 1)/2)
    ret /= sqrt(dof*pi) * gamma(dof/2)
    ret *= (1 + (x**2/dof))**(-(dof+1)/2)
    return ret

print(integrate(lambda x: t(x, 9), 1.1))