Numerically Stable Evaluation of $\frac{x^a - 1}{x^a (x - 1)}$ 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$?
 A: 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.
A: 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}$$
A: 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.
