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

Question

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

Answer 1

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.

Answer 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.

Answer 3

Here is C++ code comparing the Mathematica computation with three more methods.

1. naive
2. automatic finite differences
3. the method proposed in another answer

One important thing to note: if you make $a$ and $b$ any closer only method 2 works. In fact, method 2 works even when $a=b$ in which case it computes the derivative.

Run the code at godbolt.org.

#include <iostream>
#include <cmath>

template<typename T>
struct FD
{   
    T a, b, c;

    FD operator+(const FD& other)
    {
        return FD{a + other.a, b + other.b, c + other.c};
    }

    FD operator*(const FD& other)
    {
        return FD{a * other.a,
                  b * other.b,
                  c * other.b + a * other.c};
    }
};

template<typename T>
double FiniteDifference(T (*)(T), T a, T b)
{
    return f(FD{a, b, 1.}).c;
}

template<typename T>
T power(T x, int n)
{
    if (n == 1)
    {
        return x;
    }
    if (n % 2)
    {
        T y = power(x, n / 2);
        return x * y * y;
    }
    else
    {
        T y = power(x, n / 2);
        return y * y;
    }
}

const int n = 40;
double a = 1.9100000000000001;
double b = 1.91;

template<typename T>
T f(T x)
{
    return power(x, n);
}

//
// Mathematica gives:
// a=1910000000000000001/1000000000000000000;
// b=191/100;
// N[(a^40-b^40)/(a-b),6]
//
// 3.65058 * 10^12 
//

int main()
{
    // Obvious method
    std::cout << (pow(a, n) - pow(b, n)) / (a - b) << std::endl;
    // Automatic finite differences
    std::cout << FiniteDifference(f, a, b) << std::endl;
    // Method proposed in other answer
    std::cout << pow(a, n) * expm1(n * log1p((b - a) / a))/(b - a) << std::endl;
}
```