# Numerically stable evaluation of $x\ln(x)$

I have numerical difficulties with the function $$x\ln(x/x_0)-x+x_0$$ with $x\ge0$ and $x_0>0$ since it bears the evil evaluation of $\ln(x)$ for $x\rightarrow0$ which is the domain of interest for my numerical application. I am thus looking for a approximate function that :

• is numerically stable for $x$ in the range of (at least) $[0,x_{max}]$, with $x_{max}$ a physical property that is real and strictly positive, say between 5 and 20, usually 10.
• has controllable error (for instance by systematically improving the approximation by increasing some order, like in a series expansion)
• must be derivable.
• must not be indetermined for $x=0$.
• must contain only one minima, like the original function.

If that can be of any help, $x_0>0$ is usually of the order of 0.033.

I tested Series expansions and Padé approximants but none fill my requirements.

I think the $-x+x_0$ part is of no-importance in all of this.

I would really love some help. Thanks!

• What exactly do you mean by "numerically unstable"? I see no evidence that the expression, as is, has poor relative error as $x$ approaches zero, except for $x$ close to $x_0$, due to subtractive cancellation. Obviously the $log$ prevents you from getting a finite result when $x=0$, which you can fix in any manner suitable for your context, which has not been described. In particular, as $x$ approaches zero, the expression approaches $x_0$. – njuffa Oct 17 '15 at 13:44
• The general problem is that for $x$ close to $x_0$, the terms $-x+x_0$ and $x \ln(x/x_0)$ are of similar magnitude but of opposite sign, which means subtractive cancellation occurs. Based on the Sterbenz lemma, $-x+x_0$ is exact for $\frac{1}{2}x_0 \lt x \lt 2x$. So if we could compute $x \ln(x/x_0)$ with increased accuracy, it would help. One way of doing this is to use a fused multiply-add (FMA) operation. For example, C and C++ offer the functions fma(), fmaf() that implement this operation. Do you have FMA at your disposal? – njuffa Oct 17 '15 at 16:31
• The accuracy of $\ln(x/x_0)$ for $x$ near $x_0$ can be improved by using the C\C++ function log1p(x) which computes $\ln(1+x)$: Instead of log(x/x0), code log1p((x0 - x) / -x0). – njuffa Oct 17 '15 at 16:55

## 2 Answers

If $x > 0$ and $b^{-n} \le x < b^{-n+1}$, where $b$ can be any of $e, 2, 10,$ or otherwise depending on your inclination, since $x = b^{-n}c$ where $1 \le c < b$, we have $\ln(x) =-n\ln(b)+\ln(c)$.

Other than for $0$ itself, the expression seems well behaved as $x\rightarrow 0$. Based on the information in the question, $\epsilon \lt x_0 \lt 1$, so there is no intermediate underflow in the computation of $x/x_0$. There is also no overflow in either the division or the product since $x \le x_{max} \approx 20$. Since, as $x \rightarrow 0$, the expression $x \ln(x/x_0) -x + x_0$ approaches $x_0$, zero could easily be handled as a special case.

However, there are numerical issues for $x$ close to $x_0$. In that case, there is subtractive cancellation, as $x \ln(x/x_0) -x + x_0 = O((x-x_0)^2)$. In addition, the accuracy of $\ln(x/x_0)$ suffers as $x \rightarrow x_0$ and thus $x/x_0 \rightarrow 1$. These issues can be substantially mitigated (but not completely eliminated) as follows.

According to the Sterbenz lemma, when binary floating-point arithmetic is used, $x_0 - x$ is exact when $\frac{1}{2}x_0 \lt x \lt 2x_0$. It is therefore advantageous to compute the logarithmic term using the function log1p() for $x$ in that interval, where log1p(x) computes $\ln(1+x)$. So instead of log (x / x0) we code log1p ((x0 - x) / -x0).

Since $x_0 - x$ is exact, with trailing zero bits, we would want to compute $x \ln(x/x_0)$ to better than working precision, thus retaining additional valid bits even though the leading bits cancel during effective subtraction. An easy way of achieving this is through the use of the fused multiply-add (FMA) operation. This operation is exposed in C/C++ as the function fma(), where fma(a,b,c) computes $ab+c$ with only a single rounding at the end.

Putting the various pieces together, we arrive at the following example implementation for C/C++. Implementations for other environments should be easy to derive from this, especially if they offer ready access to fma() and log1p() functionality.

/* computes x*log(x/x0)-x+x0 for x >= 0, eps < x0 < 1 */

double func (double x, double x0)
{
double y;
if (x == 0) { // zero
y = x0;
} else if ((x <= 0.5 * x0) || (x >= 2.0 * x0)) { // x not close to x0
y = fma (x, log (x / x0), x0 - x);
} else {
y = fma (x, log1p ((x0 - x) / -x0), x0 - x); // x close to x0
}
return y;
}