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

Question

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!

Answer 1

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

Answer 2

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;
}

To achieve (close to) full machine accuracy, one suitable, if somewhat slow approach is the use of a Taylor series expansion when $x$ is very close to $x_0$. Using the command taylor series x*ln(x/k)-x+k at x=k at Wolfram Alpha gives us the following expansion around $x_0$:

$${x\log\frac{x}{x_0}-x+x_0}=\sum_{n=1}^{\infty}{\frac{(x_0-x)^{(n+1)}}{n(n+1)x_0}}; \left|x_0-x\right|\lt x_0$$

By adjusting the number of terms of the expansion used in the computation one can achieve varying levels of accuracy. In the code below I maximized accuracy by adding terms until the size of the last term falls below machine round-off. When using $\epsilon \lt x_0 \lt 1$ and $0 \le x \le 20$, the maximum errors observed with IEEE-754 double-precision computation with hundreds of millions of test vectors are on the order of 8 ulps.

#include <float.h>
#include <math.h>

double func (double x, double x0)
{
    double y;
    if (x == 0) { // x is 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 if ((x <= 0.7 * x0) || (x >= 1.4 * x0)) { // x quite close to x0
        y = fma (x, log1p ((x0 - x) / -x0), x0 - x);
    } else { // x very close to x0. Use series expansion around x0
        // sum_{n=1}^{inf} (x0-x)**(n+1) / (x0*n*(n+1)); |x0-x| < x0
        const double eps = DBL_EPSILON; // 2**(-52)
        const double z = x0 - x;
        const double f = z / x0;
        double t, q = z, n = 1;
        y = 0;
        do {
            q = q * f;
            t = q / (n * (n + 1));
            y = y + t;
            n = n + 1;
        } while (fabs (t) > (eps * fabs (y)));
    }
    return y;
}