# Method to compute $a^n - b^n$

Given two floating point numbers $$a,b$$ with $$a > b$$ and an integer $$n$$, what is the most accurate way to compute $$a^n - b^n$$ ? We can assume both $$a,b$$ are between 1 and 2. Lets assume both $$a^n$$ and $$b^n$$ can be represented in floating point arithmetic. If I compute the expression as written would I lose accuracy if $$n$$ is large? Or would it be more accurate to do something like: $$a^n - b^n = a^n \left(1 - \left( \frac{b}{a} \right)^n \right)$$

• If it helps we can assume $a,b$ are in $[1,2]$. The computation comes from computing an integral $\int_b^a x^{n-1} dx$
– vibe
Jan 14 at 4:24
• You can try to compute the integral numerically. Jan 14 at 4:28
• @lightxbulb so you mean to use a Gauss-Legendre quadrature or similar? Rather than directly evaluating the analytic solution?
– vibe
Jan 14 at 4:42

That computation is ill-conditioned anyway when $$a$$ and $$b$$ are close. This is not something that you can fix by switching to a different method: any method that uses floating-point computations will be safe to use only if you have exact values of the inputs $$a$$ and $$b$$ available as floating-point numbers. So even working with $$a=\frac43$$ is dangerous, because it is not an exactly representable floating point value.
If you can, you should try to avoid this computation altogether and reframe your computation so that you only use $$\frac{a^n-b^n}{a-b},$$ which can be computed accurately (at least when $$a,b \in [1,2]$$) as $$a^{n-1} + a^{n-2}b + \dots + ab^{n-2} + b^{n-1}$$ because all the summands have the same sign and there is no cancellation.
you got pretty close to a good method. Specifically what you want to do is compute a^n*exp(n*log1p(-b/a)) when the values are close together. this will avoid catastrophic cancellation.
• This is wrong. It should be -a^n*expm1(n*log1p((b-a)/a)). In the example this gives 7.526705395707148e+24, which is close to the exact result in the limitations of floating-point numbers. Jan 14 at 7:06