About faster approximation of log(x)

Question

I had written a code a while ago which attempted to calculate $log(x)$ without using library functions. Yesterday, I was reviewing the old code, and I tried to make it as fast as possible, (and correct). Here's my attempt so far:

const double ee = exp(1);

double series_ln_taylor(double n){ /* n = e^a * b, where a is an non-negative integer */
    double lgVal = 0, term, now;
    int i, flag = 1;

    if ( n <= 0 ) return 1e-300;
    if ( n * ee < 1 )
        n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */

    for ( term = 1; term < n ; term *= ee, lgVal++ );
    n /= term;

    /* log(1 - x) = -x - x**2/2 - x**3/3... */
    n = 1 - n;
    now = term = n;
    for ( i = 1 ; ; ){
        lgVal -= now;
        term *= n;
        now = term / ++i;
        if ( now < 1e-17 ) break;
    }

    if ( flag == -1 ) lgVal = -lgVal;

    return lgVal;
}

Here I am trying to find $a$ so that $e^{a}$ is just over n, and then I add the logarithm value of $\frac{n}{e^{a}}$, which is less than 1. At this point, the Taylor expansion of $log(1\ - \ x)$ can be used without worrying.

I have recently grown an interest in numerical analysis, and that's why I can't help asking the question, how much faster this code segment can be run in practice, while being correct enough? Do I need to switch to some other methods , for example, using continued fraction, like this?

The $log(x)$ function provided with C standard library is almost 5.1 times faster than this implementation.

UPDATE 1: Using the hyperbolic arctan series mentioned in Wikipedia, the computation seems to be almost 2.2 times slower than C standard library log function. Though, I have not extensively checked the performance, and for larger numbers, my current implementation seems to be REALLY slow. I want to check both of my implementation for error bound and average time for a wide range of numbers if I can manage. Here is my second effort.

double series_ln_arctanh(double n){ /* n = e^a * b, where a is an non-negative integer */
    double lgVal = 0, term, now, sm;
    int i, flag = 1;

    if ( n <= 0 ) return 1e-300;
    if ( n * ee < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */

    for ( term = 1; term < n ; term *= ee, lgVal++ );
    n /= term;

    /* log(x) = 2 arctanh((x-1)/(x+1)) */
    n = (1 - n)/(n + 1);

    now = term = n;
    n *= n;
    sm = 0;
    for ( i = 3 ; ; i += 2 ){
        sm += now;
        term *= n;
        now = term / i;
       if ( now < 1e-17 ) break;
    }

    lgVal -= 2*sm;

    if ( flag == -1 ) lgVal = -lgVal;
    return lgVal;
}

Any suggestion or criticism is appreciated.

UPDATE 2: Based on suggestions made below, I have added some incremental changes here, which is around 2.5 times slower than standard library implementation. However, I have tested it only for integers $\leq1e^{8}$ this time, for larger numbers the runtime would increase. For now. I do not yet know of techniques to generate random double numbers $\leq1e^{308}$, so it's not yet fully benchmarked. To make the code more robust, I have added corrections for corner cases. The average error for the tests I made is around $4e^{-15}$.

double series_ln_better(double n){ /* n = e^a * b, where a is an non-negative integer */
    double lgVal = 0, term, now, sm;
    int i, flag = 1;

    if ( n == 0 ) return -1./0.; /* -inf */
    if ( n < 0 ) return 0./0.;   /* NaN*/
    if ( n < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */

    /* the cutoff iteration is 650, as over e**650, term multiplication would
       overflow. For larger numbers, the loop dominates the arctanh approximation
       loop (with having 13-15 iterations on average for tested numbers so far */

    for ( term = 1; term < n && lgVal < 650 ; term *= ee, lgVal++ );
    if ( lgVal == 650 ){
        n /= term;
        for ( term = 1 ; term < n ; term *= ee, lgVal++ );
    }
    n /= term;

    /* log(x) = 2 arctanh((x-1)/(x+1)) */
    n = (1 - n)/(n + 1);

    now = term = n;
    n *= n;
    sm = 0;

    /* limiting the iteration for worst case scenario, maximum 24 iteration */
    for ( i = 3 ; i < 50 ; i += 2 ){
        sm += now;
        term *= n;
        now = term / i;
        if ( now < 1e-17 ) break;
    }

    lgVal -= 2*sm;

    if ( flag == -1 ) lgVal = -lgVal;

    return lgVal;
}

Answer 1

This is not really an authoritative answer, more a list of issues I think you should consider, and I haven't tested your code.

0. How did you test your code for correctness and speed? Both are important, somewhat tricky to do well, and you don't give details. In other words, if I compare your function with the log on my machine, will I also get the same numbers, $2.1$, $5.1$? In my experience with reading other people's timing benchmarks in academic literature, it requires a lot of care and precision to get reproducible timings, which are the only timings anyone is ever going to care about. Microbenchmarks especially are notoriously unreliable.

1. The common problem with evaluating a function $f(x)$ directly with its (unmodified) Taylor series is the number of terms necessary for convergence. There are 52 bits in the mantissa of a double, so when $n\approx\frac12$ at the beginning of the Taylor series loop, you can expect the loop to take about 50-ish iterations. That's pretty expensive and should be optimized.

1.5. Have you checked your code for large $n$? I tried 1.7976e+308, which leads to term=inf, and then to n=1 in the Taylor series loop, leading to extremely slow convergence: it converges like the Harmonic series, i.e. it doesn't but there will be at most $10^{17}$ terms. As a rule of thumb, you should have some kind of a "max-iteration-count" bound for the loop. In this case it behaves like this because $n$ is finite, but term *= e overflows to $\infty$ in the argument reduction loop. The correct answer is $709.78266108405500745\ldots$.

2. Functions implemented in standard libraries are expected to be extremely robust. Returning $10^{-300}$ ($\approx0$) as the logarithm of a negative (or zero) number is not correct. The logarithm of $0$ should be $-\infty$, logarithm of a negative number should be NaN.

I suspect that with a little effort you can sacrifice some of that robustness for performance, e.g., by restricting the argument range or returning slightly less accurate results.

3. Performance of this kind of code can depend a lot on the CPU architecture it's running on. It's a deep and involved topic, but CPU manufacturers like Intel publish optimization guides that explain the different interactions between your code and the CPU it's running on. Caching can be relatively straightforward, but things like branch prediction, instruction-level parallelism, and pipeline stalls due to data dependencies are difficult to see precisely in high-level code, but matter a lot for performance.

4. Implementing a function like this correctly usually means that you guarantee that for the input floating-point number $\tilde x$, the output $\tilde y = \tilde f(\tilde x)$ is within a certain distance of the nearest floating-point number to the true value $y = f(\tilde x)$. Verifying this is not totally trivial, there isn't any evidence in your code that you've done this, so I don't know if your function is correct (I'm sure it's quite accurate, but how accurate?). This is not the same as showing that the Taylor series converges, due to the presence of floating-point roundoff errors.

4.5. A good way to test an untested function for accuracy, would be to evaluate it at each one of the four billion (fewer if you are doing argument reduction correctly, as here) single-precision floats, and compare the errors with the standard log from libm. Takes a bit of time, but at least it's thorough.

5. Because you know from the start the precision of doubles, you don't have to have an unbounded loop: the number of iterations can be figured out up front (it's probably about 50). Use this to remove branches from your code, or at least set the number of iterations in advance.

All the usual ideas about loop unrolling apply too.

6. It is possible to use approximation techniques other than Taylor series. There are also Chebyshev series (with the Clenshaw recurrence), Pade approximants, and sometimes root-finding methods like Newton's method whenever your function can be recast as the root of a simpler function (e.g., the famous sqrt trick).

Continued fractions are probably not going to be too great, because they involve division, which is much more expensive than multiplies/adds. If you look at _mm_div_ss at https://software.intel.com/sites/landingpage/IntrinsicsGuide/, division has latency of 13-14 cycles and throughput of 5-14, depending on architecture, compared with 3-5/0.5-1 for multiply/add/madd. So in general (not always) it makes sense to try to eliminate divisions as much as possible.

Unfortunately, mathematics isn't such a great guide here, because expressions with short formulas are not necessarily the fastest ones. Mathematics doesn't penalize divisions, for example.

7. Floating-point numbers are stored internally in the form $x = m\times 2^e$ (mantissa $m$, $\frac12<m\leq1$, exponent $e$). The natural log of $x$ is much less natural than the base-2 log, for which the first part of your code can be replaced with one call to frexp.

8. Compare your log with the log in libm or openlibm(e.g.: https://github.com/JuliaLang/openlibm/blob/master/src/e_log.c). This is by far the easiest way to find out what other people have already figured out. There are also specially-optimized versions of libm specific to CPU manufacturers, but those usually don't have their source code published.

Boost::sf has some special functions, but not the basic ones. It might be instructive to look at the source of log1p, though: http://www.boost.org/doc/libs/1_58_0/libs/math/doc/html/math_toolkit/powers/log1p.html

There are also open-source arbitrary-precision arithmetic libraries like mpfr, which might use different algorithms than libm due to the higher precision required.

9. Higham's Accuracy and Stability of Numerical Algorithms is a good upper-level introduction to analyzing errors of numerical algorithms. For approximation algorithms themselves, Approximation Theory Approximation Practice by Trefethen is a good reference.

10. I know this is said a bit too often, but reasonably large software projects rarely depend on runtime of one small function being called over and over again. It's not so common to have to worry about performance of log, unless you've profiled your program and made sure it's important.

Answer 2

Kirill's answer already touched on a large number of relevant issues. I would like to expand on some of them based on practical math library design experience. A note up front: math library designers tend to use every published algorithmic optimization, as well as many machine specific optimizations, not all of which will be published. The code frequently will be written in assembly language, rather than use compiled code. It is therefore unlikely that a straightforward and compiled implementation will achieve more than 75% of the performance of an existing high-quality math library implementation, assuming identical feature sets (accuracy, special cases handling, error reporting, rounding-mode support).

In terms of accuracy, for elementary functions (e.g. $exp$, $log$) as well as simple special functions (e.g. $erfc$, $\Gamma$), ulp error has replaced relative error as the relevant error metric. A common design goal for elementary functions is a maximum error of less than 1 ulp, resulting in a faithfully-rounded function. A faithfully rounded function returns either the floating-point number closest to the infinitely precise result, or an immediately adjacent floating-point number.

Accuracy is typically assessed by comparison with a (third-party) higher-precision reference. Single-argument single-precision functions can easily be tested exhaustively, other functions require testing with (directed) random test vectors. Clearly one cannot compute infinitely precise reference results, but research into the Table-Maker's Dilemma suggests that for many simple functions it suffices to compute a reference with a precision of about three times the target precision. See for example:

Vincent Lefèvre, Jean-Michel Muller, "Worst Cases for Correct Rounding of the Elementary Functions in Double Precision". In Proceedings 15th IEEE Symposium on Computer Arithmetic, 2001,111-118).(preprint online)

In terms of performance, one has to distinguish between optimizing for latency (important when one looks at the execution time of dependent operations), versus optimizing for throughput (relevant when considering execution time of independent operations). Over the past twenty years, the proliferation of hardware parallelization techniques such as instruction level parallelism (e.g. superscalar, out-of-order processors), data-level parallelism (e.g. SIMD instructions), and thread-level parallelism (e.g. hyper-threading, multi-core processors) has lead to an emphasis on computational throughput as the more relevant metric.

Core approximations for elementary functions are almost exclusively minimax approximations. Due to the relatively high cost of division, polynomial minimax approximations. Tools such as Mathematica or Maple have built-in facilities for generating these; there are also specialized tools such as Sollya. For the logarithm, the basic core approximation choices, after argument reduction to values close to unity, are $log (1+x) = p(x)$ and $log (x) = 2 \cdot atanh ((x-1)/ (x+1)) = p(((x-1)/ (x+1))^{2})$, where $p$ is a polynomial minimax approximation. Performance-wise, the former is usually preferred for single-precision implementations (see this answer for a worked example), while the latter is preferred for double-precision implementations.

The fused multiply-add operation (FMA), first introduced by IBM 25 years ago, and now available on all major processor architectures, is a crucial building block of modern math library implementations. It provides rounding error reduction, gives limited protection against subtractive cancellation, and vastly simplifies double-double arithmetic.

The exemplary IEEE-754 double-precision C99 implementation of log() below demonstrates the use of FMA (exposed in C99 as the fma() standard math function), along with very limited use of double-double techniques to improve the accuracy of products with transcendental constants. Two different core approximations are provided, both delivering faithfully-rounded results, as demonstrated by testing with $2^{33}$ random test vectors. The minimax approximations used were computed with my own tools based on the Remez algorithm.

#include <math.h>

/* compute natural logarithm

   USE_ATANH == 1: maximum error found: 0.83482 ulp @ 0.7012829191167614
   USE_ATANH == 0: maximum error found: 0.83839 ulp @ 1.2788954397331760
*/
double my_log (double a)
{
    const double LOG2_HI = 0x1.62e42fefa39efp-01; // 6.9314718055994529e-01
    const double LOG2_LO = 0x1.abc9e3b39803fp-56; // 2.3190468138462996e-17
    double m, r, i, s, t, p, f, q;
    int e;

    m = frexp (a, &e);
    if (m < 0.70703125) { // 181/256
        m = m + m;
        e = e - 1;
    }
    i = (double)e;

    /* m in [181/256, 362/256] */

#if USE_ATANH
    /* Compute q = (m-1) / (m+1) */
    p = m + 1.0;
    m = m - 1.0;
    q = m / p;

    /* Compute (2*atanh(q)/q-2*q) as p(q**2), q in [-75/437, 53/309] */
    s = q * q;
    r =             0x1.2f1da230fb057p-3;  // 1.4800574027992994e-1
    r = fma (r, s,  0x1.399f73f934c01p-3); // 1.5313616375223663e-1
    r = fma (r, s,  0x1.7466542530accp-3); // 1.8183580149169243e-1
    r = fma (r, s,  0x1.c71c51a8bf129p-3); // 2.2222198291991305e-1
    r = fma (r, s,  0x1.249249425f140p-2); // 2.8571428744887228e-1
    r = fma (r, s,  0x1.999999997f6abp-2); // 3.9999999999404662e-1
    r = fma (r, s,  0x1.5555555555593p-1); // 6.6666666666667351e-1
    r = r * s;

    /* log(a) = 2*atanh(q) + i*log(2) = LOG2_LO*i + p(q**2)*q + 2q + LOG2_HI*i.
       Use K.C. Ng's trick to improve the accuracy of the computation, like so:
       p(q**2)*q + 2q = p(q**2)*q + q*t - t + m, where t = m**2/2.
    */
    t = m * m * 0.5;
    r = fma (q, t, fma (q, r, LOG2_LO * i)) - t + m;
    r = fma (LOG2_HI, i, r);

#else // USE_ATANH

    /* Compute f = m -1 */
    f = m - 1.0;
    s = f * f;

    /* Approximate log1p (f), f in [-75/256, 106/256] */
    r = fma (-0x1.961d64ddd82b6p-6, f, 0x1.d35fd598b1362p-5); // -2.4787281515616676e-2, 5.7052533321928292e-2
    t = fma (-0x1.fcf5138885121p-5, f, 0x1.b97114751d726p-5); // -6.2128580237329929e-2, 5.3886928516403906e-2
    r = fma (r, s, t);
    r = fma (r, f, -0x1.b5b505410388dp-5); // -5.3431043874398211e-2
    r = fma (r, f,  0x1.dd660c0bd22dap-5); //  5.8276198890387668e-2
    r = fma (r, f, -0x1.00bda5ecdad6fp-4); // -6.2680862565391612e-2
    r = fma (r, f,  0x1.1159b2e3bd0dap-4); //  6.6735934054864471e-2
    r = fma (r, f, -0x1.2489f14dd8883p-4); // -7.1420614809115476e-2
    r = fma (r, f,  0x1.3b0ee248a0ccfp-4); //  7.6918491287915489e-2
    r = fma (r, f, -0x1.55557d3b497c3p-4); // -8.3333481965921982e-2
    r = fma (r, f,  0x1.745d4666f7f48p-4); //  9.0909266480136641e-2
    r = fma (r, f, -0x1.999999d959743p-4); // -1.0000000092767629e-1
    r = fma (r, f,  0x1.c71c70bbce7c2p-4); //  1.1111110722131826e-1
    r = fma (r, f, -0x1.fffffffa61619p-4); // -1.2499999991822398e-1
    r = fma (r, f,  0x1.249249262c6cdp-3); //  1.4285714290377030e-1
    r = fma (r, f, -0x1.555555555f03cp-3); // -1.6666666666776730e-1
    r = fma (r, f,  0x1.999999999759ep-3); //  1.9999999999974433e-1
    r = fma (r, f, -0x1.fffffffffff53p-3); // -2.4999999999999520e-1
    r = fma (r, f,  0x1.555555555555dp-2); //  3.3333333333333376e-1
    r = fma (r, f, -0x1.0000000000000p-1); // -5.0000000000000000e-1

    /* log(a) = log1p (f) + i * log(2) */
    p = fma ( LOG2_HI, i, f);
    t = fma (-LOG2_HI, i, p);
    f = fma ( LOG2_LO, i, f - t);
    r = fma (r, s, f);
    r = r + p;
#endif // USE_ATANH

    /* Handle special cases */
    if (!((a > 0.0) && (a <= 0x1.fffffffffffffp1023))) {
        r = a + a;  // handle inputs of NaN, +Inf
        if (a  < 0.0) r =  0.0 / 0.0; //  NaN
        if (a == 0.0) r = -1.0 / 0.0; // -Inf
    }
    return r;
}