Numerically Stable Evaluation of $\frac{x^a - 1}{x^a (x - 1)}$

Question

I want to evaluate the following function as part of a computer program: $$f(x) = \frac{x^a - 1}{x^a (x - 1)}, x \in \mathbb{R}^+, a \in \mathbb{R}_0^+$$

Mathematically, this function is fairly well-behaved, with the singularity at $x = 1$ being a removable singularity. However, even though the singularity can be treated separately, the evaluation is numerically unstable in the neighborhood of $x = 1$ due to both the numerator and the denominator becoming very small.

Even though $x = 1$ is a root of the numerator, I have been unable to algebraically express the numerator as a product $(x - 1)r$ in order to cancel the linear factor $x - 1$ to remove the singularity. I know it is possible for $a \in \mathbb{N}$ via the formulas for the geometric progression, but the resulting polynomial is expensive to evaluate and in my case, $a$ can be an arbitrary nonnegative real number, so this does not apply anyways.

Is there a way to algebraically express a function that is $f$ but with the singularity removed? If not, is there a way to evaluate $f$ while avoiding the numerical instability around $x = 1$?

Answer 1

I want to highlight Lutz Lehmann's comment on the expm1 function which is ubiquitous. It's purpose is to efficiently compute $e^x -1$ for $x \approx 0$.

Lutz gives the formula

expm1(a*log(x))/pow(x,a)/(x-1)

which looks correct to me, and does indeed appear to be stable for $x \approx 1$. You must check for $x = 1$ to avoid division by zero.

I suspect that all the answers using Taylor series are essentially reinventing the wheel here.

Answer 2

For convenience set $x=1+s$, then near $s=0$ you can use the Taylor series

$$ \eqalign{f(1+s) &= \sum_{n=0}^\infty \frac{(-1)^n \Gamma(n+a+1)}{\Gamma(a) (n+1)!} s^n \cr &= a - \dfrac{a(a+1)}{2} s + \dfrac{a(a+1)(a+2)}{6} s^2 - \dfrac{a(a+1)(a+2)(a+3)}{24} s^3 + \ldots}$$

Answer 3

Presumably, the problems are not caused because the values near $1$ are small or large pre se, but rather by cancellation when $x\approx 1$.$\def\e{\varepsilon}$

So when $x\approx 1$, you can try to approximate $f$ by a simpler function like Taylor expansion around $1$, for example

$$f(x)\approx \frac ax-\frac12 a(a-1)(x-1)$$

When $x\approx 1$ let's write it as $x=1+\e$:

$$\begin{align} f(x)=\frac{x^a-1}{x^a(x-1)} &= \frac1{x^a}\underbrace{\frac{(1+\e)^a-1^a}{\e}}_{\textstyle\approx g'(1)+\dfrac\e2 g''(1)} \\ \end{align}$$ where $g(x) = x^a$ and thus $$g'(x) = ax^{a-1}$$ $$g''(x) = a(a-1)x^{a-2}$$

Here is a Desmos plot zoomed at $x\approx 1$. Notice that the functions in the plot have actually $a$ subtracted so that they all hit $(1,0)$ irrespective of the value of $a$. You can play around with $a$ by dragging the slider.

The red function is $f(x)$ that goes bonkers close to 1. The violet function is the approximation near 1.

As far as I know, Desmos is using IEEE-754 double, i.e. "ordinary" 64-bit binary floating point, so you get some visual quality estimates.

Answer 4

Using a simple Padé approximant

$$\frac{x^a - 1}{x^a (x - 1)}=a \frac{1+\frac{7-a}{10} (x-1)+\frac{(a-1)(a-2)}{60} (x-1)^2 } {1+\frac{2(a+3)}{5} (x-1)+\frac{(a+2)(a+3)}{20} (x-1)^2 }$$

Trying for $a=\pi$ and $x=1+10^{-k}$, some results

$$\left( \begin{array}{ccc} k & \text{approximation} & \text{exact value} \\ 1 & \color{red}{2.58756}3328300637403099800 & 2.587562502050549708509754 \\ 2 & \color{red}{3.0776347615}76885361735203 & 3.077634761563679134659739 \\ 3 & \color{red}{3.13509818768015}6052084880 & 3.135098187680155913340110 \\ 4 & \color{red}{3.14094220521706779548}2107 & 3.140942205217067795480713 \\ 5 &\color{red}{ 3.141527598719473976878852} & 3.141527598719473976878852 \end{array} \right)$$ and we could do much better.