Question: Given a real function $f$, how can one efficiently evaluate $\int\_0^{a\_i}f(x)dx$ for millions of different $a\_i$?

Applying a standard [quadrature](https://en.wikipedia.org/wiki/Numerical_integration) method (such as [Simpson's rule](https://en.wikipedia.org/wiki/Simpson%27s_rule) or [Gaussian quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature)) will incur an independent cost for each $a\_i$, since in general $f(x)$ will need to be evaluated at an entirely different set of points for each $a\_i$.

I would prefer to pre-generate/calculate a data structure at substantial initial cost, if the integral can then be evaluated from this data structure cheaply, for different $a\_i$. 

One possibility I see is to evaluate the integral at a regular array of points $x\_j=a\_{max}j/n$ (which would allow reuse of points in Simpson's method), and then to "correct" this integral on final evaluation with a small $\int\_{a\_{max}j/n}^{a\_i}f(x)dx$ which can be evaluated with a very low-order Gaussian quadrature because $j$ can be chosen to make this correction interval small. Although probably accurate, it seems like this method may give rise to discontinuities at the half-way points between subsequent $j$.

Another possibility may be to fit cubic splines to a regular array of points (evaluated from Simpson's method) and then rely on these to interpolate the integral.

I would rather follow an established approach if there is one.

Answer: 
Let's say you've ordered the points so that $0<a\_1<a\_2<\ldots<a\_n$. Then you have of course that
$$
 \int\_0^{a\_2} f(x) dx
 =
 \int\_0^{a\_1} f(x) dx +
 \int\_{a\_1}^{a\_2} f(x) dx.
$$
In other words, if you have previously computed the integrals from zero to $a\_1$ through some numerical method, then computing the integral from zero to $a\_2$ can reuse this information, and all you have to do is now compute the integral from $a\_1$ to $a\_2$, which you can do independently.

The same applied for all of the other integrals, of course.
The standard method is [double-exponential sinc quadrature](http://www.ams.org/journals/mcom/2005-74-250/S0025-5718-04-01724-7/S0025-5718-04-01724-7.pdf). The fast convergence requires some regularity assumptions on the integrand, namely that it is the restriction to the real line of a complex differentiable function bounded on a Bernstein ellipse.