Building a ln function from scratch

Question

I need to write a natural logarithm with just limited available operations. This is to be implemented as a function into a C program. The problem is, our product is security relevant so the source code has to be well documented and tested. The testing process is why we are writing all code from scratch, as the C standard librarys are not well documented and writing our own tests for it that cover the whole library would take longer as just developing what we need our self.

Long story short:

I reviewed some opensource math librarys and just encoutnered for every ln()-function/method a quite big source with some nested dependencies to other even biger math-functions.

double log(double x)
{
    int e;
    double y, z;
    /* Test for domain */
    if( x <= 0.0 )
        return 0;

    /* separate mantissa from exponent */
    /* Note, frexp is used so that denormal numbers
     * will be handled properly.
     */
    x = frexp( x, &e );
    /* logarithm using log(x) = z + z**3 P(z)/Q(z),
     * where z = 2(x-1)/x+1)
     */

    if( (e > 2) || (e < -2) )
    {
        if( x < SQRTH )
        { /* 2( 2x-1 )/( 2x+1 ) */
            e -= 1;
            z = x - 0.5;
            y = 0.5 * z + 0.5;
        }   
        else
        { /*  2 (x-1)/(x+1)   */
            z = x - 0.5;
            z -= 0.5;
            y = 0.5 * x  + 0.5;
        }

        x = z / y;
        /* rational form */
        z = x*x;
        z = x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );
        y = e;
        z = z - y * 2.121944400546905827679e-4;
        z = z + x;
        z = z + e * 0.693359375;
        goto ldone;
    }
    /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */

    if( x < SQRTH )
    {
        e -= 1;
        x = ldexp( x, 1 ) - 1.0; /*  2x - 1  */
    }   
    else
    {
        x = x - 1.0;
    }
    /* rational form */
    z = x*x;
    y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 5 ) );

    if( e )
        y = y - e * 2.121944400546905827679e-4;
    y = y - ldexp( z, -1 );   /*  y - 0.5 * z  */
    z = x + y;
    if( e )
        z = z + e * 0.693359375;

ldone:

    return( z );
}

This is code was allready hardened by me. But it is still depending on 4 other math functions (ldexp, polevl, p1evl and frexp) which I would have to implement my self aswell.

My math knowledge is allready overstrained with this code, so I couldn't verify even this is working correctly after I trimmed the code not mentioning implementing the additional implementations or the moment I had to document and jsutify the tests.

So why I'm asking here and not on StackOverflow:

I want to know, are there any easyer implementations or approximation procedures to determine the ln in a not that complex way, so I don't have to waste alot of resources for studying the maths behind?

Answer 1

For functions such as $\log$ from a standard library, I think it would be a mistake to reimplement them from scratch in the hope that your own implementation would be somehow "better". The standard implementations have been optimized for speed and accuracy for several decades, and it's unlikely that you can improve on it unless you are a highly specialized expert both in numerical analysis and the relevant hardware architecture.

That said, there are other libraries that are in some sense better, for example returning values that are guaranteed to be correctly rounded. For example, you might be interested in crlibm

---

Added: If you you want something simpler, and can live with a limited range of allowed input values, you can for example find a rational approximation of $\log x$. For example, computing the "(3,3) Padé approximation" of $\log x$ around $x=1$ gives: $$ \log x \approx \frac{\frac{11}{60}(x-1)^3 + (x-1)^2 + (x-1)}{-\frac12 + \frac32 x + \frac 35(x-1)^2 + \frac1{20}(x-1)^3} $$ giving reasonable accuracy, at least on the interval $[0.5, 4]$. Here is a graph of the two functions

and a graph of the quotient (showing relative accuracy)

Answer 2

I had to implement logarithms on rational numbers.

This is how I did it:

Occording to Wikipedia, there is the Halley-Newton approximation method

which can be used.

Using Newton's method, the iteration simplifies to (implementation), which has cubic convergence to ln(x).

// Using Newton's method, the iteration simplifies to (implementation) 
// which has cubic convergence to ln(x).
public static double ln(double x, double epsilon)
{
    double yn = x - 1.0d; // using the first term of the taylor series as initial-value
    double yn1 = yn;
do
{
    yn = yn1;
    yn1 = yn + 2 * (x - System.Math.Exp(yn)) / (x + System.Math.Exp(yn));
} while (System.Math.Abs(yn - yn1) > epsilon);

return yn1;

}

This is not C, but C#, but I'm sure anybody capable to program in C will be able to deduce the C-Code from that.

Forthermore, since

logn(x) = ln(x)/ln(n).

You have therefore just implemented logN as well.

public static double log(double x, double n, double epsilon)
{
    return ln(x, epsilon) / ln(n, epsilon);
}

where epsilon (error) is the minimum precision.

Now as to speed, you're probably better of using the ln-cast-in-hardware, but as I said, I used this as a base to implement logarithms on a rational numbers class.

Arbitrary precision might be more important than speed, under certain circumstances.

Then, use the logarithmic identities for rational numbers:
logB(x/y) = logB(x) - logB(y)