6
$\begingroup$

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) $$

$\endgroup$
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$
    – vibe
    Jan 14, 2023 at 4:24
  • 2
    $\begingroup$ You can try to compute the integral numerically. $\endgroup$
    – lightxbulb
    Jan 14, 2023 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$
    – vibe
    Jan 14, 2023 at 4:42

2 Answers 2

6
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
2
  • 5
    $\begingroup$ 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. $\endgroup$ Jan 14, 2023 at 7:06
  • 3
    $\begingroup$ that will teach me to not answer these at 2am on a phone :) $\endgroup$ Jan 14, 2023 at 7:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.