# Efficient way of calculating a cumulative integral with prefactor

I have a grid of points $$x_i$$ and corresponding function values $$y_i=f(x_i)$$. I'm interesting in something like the cumulant of $$f$$, but it has an awkward prefactor. The desired quantity we'll call

$$z(x) = e^{x}\int_{x}^{\infty}e^{-q}f(q)dq$$

And this quantity is needed for each grid point, $$z_i=z(x_i)$$

Let us assume for numerical sake that $$\infty$$ is really a sufficiently large $$q$$ such that $$f(q)$$ goes smoothly to zero, and that $$f(q)$$ is a nicely behaved, non-negative function for any $$q$$ we might care about. In this particular context, its a probability distribution. Note that I don't always have a functional form for $$f$$, so an analytical solution isn't particularly helpful here.

Were it not for the $$e^x$$ prefactor, integrating this would be rather simple, numerically. Either use scipy.integrate.cumulative_trapezoid or do the trapezoidal rule yourself by taking

integrand = np.exp(-x) * y
almostZ = dx * (np.cumsum(integrand[::-1])[::-1] - integrand/2 - integrand[-1]/2)
Z = almostZ * np.exp(x)


where dx is the grid spacing for the $$x_i$$. Done this way, you get $$z_i$$ in roughly $$O(N)$$ scaling, which is great.

The problem arises when $$x$$ gets sufficiently large that either $$e^x$$ or $$e^{-x}$$ are up against the limits of numerical precision, even though the product of the two terms is a reasonably sized number. Doing this in a loop is the current state of the "answer":

for i in range(len(x)):
Z[i] = np.trapz(np.exp(x[i] - x[i:]) * y[i:], x[i:])


This ensures that there's no really huge numbers multiplied by really small numbers, and is instead a sum of reasonably sized numbers. But it also involves a lot of double counting and $$O(N^2)$$-ness, which is a bit of a pain when this is something that's being calculated as integrating an ODE (via solve_ivp).

Is there a trick to doing this sort of computation, or am I just doomed to $$O(N^2)$$?

• Would arbitrary precision arithmetic help? The individual operations would likely be slower, but it would allow you to use your proposed O(N) algorithm. I suppose it would depend how large a system you are considering whether this would be worth it. Jan 11 at 2:32
• @Tyberius I suppose that's the obvious answer. I totally didn't think of that, and lo and behold it actually helps quite a bit. The overhead for the arbitrary precision hurts a little at small system sizes but at the number of points I need for $x_i$ its a 30% increase in speed! Jan 12 at 21:36