from math import *

# Define the Gauss[7] nodes and weights
gauss_nodes = [ \
 -0.949107912342759, -0.741531185599394, -0.405845151377397, 0.0 ,\
  0.405845151377397,  0.741531185599394,  0.949107912342759]

gauss_weights = [ \
 0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, \
 0.381830050505119, 0.279705391489277, 0.129484966168870]

# Define the Kronrod[15] nodes and weights
kronrod_nodes = [ \
 -0.991455371120813, -0.949107912342759, -0.864864423359769, -0.741531185599394,\
 -0.586087235467691, -0.405845151377397, -0.207784955007898,  0.0, \
  0.207784955007898,  0.405845151377397,  0.586087235467691,  0.741531185599394, \
  0.864864423359769,  0.949107912342759,  0.991455371120813]

kronrod_weights = [ \
 0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, \
 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, \
 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525, \
 0.104790010322250, 0.063092092629979, 0.022935322010529]


# Define the Gauss-Kronrod[7,15] quadrature
def gauss_kronrod_7_15(f):
    "Gauss-Kronrod 7-15 Quadrature"
    gauss, kronrod = 0., 0.
    for n,w in zip(gauss_nodes, gauss_weights):
        gauss += f(n)*w

    for n,w in zip(kronrod_nodes, kronrod_weights):
        kronrod += f(n)*w
		 
    err = (200*abs(gauss - kronrod))**(1.5)
    return (kronrod, err)


# Define the standard normal distribution N(0,1)
def normalDist(x):
    "N(0,1)"
    return 1./sqrt(2.*pi)*exp(-(x**2)/2.)


# Test the quadrature
print " GK[7,15] of N(0,1) is (value, err):",
print gauss_kronrod_7_15(normalDist)
print "\n (Maple gives 0.6826894921370859)"