# Integration of an exponential function over a triangle. Numerically stable approximation to the infinite sum?

This question is motivated by my need to evaluate the following integral for general a and b, \begin{align} I(a, b) &= \int_{0}^{1}\mathrm{d}x\;\int_{0}^{1-x}\mathrm{d}y \left( e^{ax + by} \right) \\ &= \frac{a (1 - e^{b}) - b(1 - e^a)}{ab(a-b)}. \end{align}

I could do cases for small a, b, and (a-b). But, I'm wondering if there is a better way. The $a$ and $b$ values are variable in a particular way,
$$a=-\sigma \alpha$$ and, $$b=-\sigma \beta.$$

The values of $\alpha$ and $\beta$ are fixed, but $\sigma$ can be any non-negative number. I want a fast rule for evaluating the integral for any new $\sigma$.

One thing that I tried was expanding the exponential functions in power series. This gives, $$I(a, b) = \frac{1}{2} + \frac{1}{a-b}\sum_{n=2}^{\infty}\frac{a^{n}-b^{n}}{(n+1)!}$$

Maybe I can cancel out the $(a-b)$ in the denominator and truncate take a partial sum? It would be great if I could get an expression like the following, $$I(-\sigma \alpha, -\sigma\beta) = c_{0}(\alpha, \beta) + c_{1}(\alpha, \beta) \sigma + c_{2}(\alpha, \beta)\sigma^{2} + \cdots,$$

so that I can just store the constants ($c_{0}$, $c_{1}$, $c_{2}$, ...) and evaluate a low order polynomial in $\sigma$.

Does anybody have any insight into how to calculate this integral quickly for new values of sigma?