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) $$
3
-
1$\begingroup$ 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$ $\endgroup$– vibeJan 14 at 4:24
-
2$\begingroup$ You can try to compute the integral numerically. $\endgroup$– lightxbulbJan 14 at 4:28
-
1$\begingroup$ @lightxbulb so you mean to use a Gauss-Legendre quadrature or similar? Rather than directly evaluating the analytic solution? $\endgroup$– vibeJan 14 at 4:42
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
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.
answered Jan 14 at 11:42
$\begingroup$
$\endgroup$
2
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.
-
5$\begingroup$ This is wrong. It should be
-a^n*expm1(n*log1p((b-a)/a)). In the example this gives7.526705395707148e+24, which is close to the exact result in the limitations of floating-point numbers. $\endgroup$ Jan 14 at 7:06 -
3$\begingroup$ that will teach me to not answer these at 2am on a phone :) $\endgroup$ Jan 14 at 7:09