Background

To approximate an integral over the reference interval [-1, 1], one arranges points x_i and weights w_i such that the number

sum(f(x_i) * w_i)

approximates the integral of some functions f (typically a set of polynomials) reasonably well. Sophisticated schemes have many points and weights, and naturally the more weights there are, the smaller each individual weight must be.

Bottom line: You'll have to sum up many small values.

When numpy.sum doesn't cut it

I noticed that numpy.sum produces large enough round-off errors for the tests to fail with as few as 7 summands (being -44.34121805, -15.72356145, -0.04563372, 0., 0.04563372, 15.72356145, 44.34121805, not shown in full precision). math.fsum always worked fine, but cannot be used anymore since I'd like to compute the sum for many intervals at once.

Work so far

Kahan's summation

def kahan_sum(a, axis=0):
    s = numpy.zeros(a.shape[:axis] + a.shape[axis+1:])
    c = numpy.zeros(s.shape)
    for i in range(a.shape[axis]):
        # http://stackoverflow.com/a/42817610/353337
        y = a[(slice(None),) * axis + (i,)] - c
        t = s + y
        c = (t - s) - y
        s = t.copy()
    return s

does the trick. Note, however, the use of a Python loop here, which leads to a performance degradation if the number of points and weights is large.

math.fsum uses a different algorithm which may or may not be better suited than Kahan.

See here for a previous discussion of the topic.