0
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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 enter image description here

and a graph of the quotient (showing relative accuracy)

enter image description here

$\endgroup$
  • $\begingroup$ I don't want to improve them, I can even live with a decent fallback in accuracee and performance isn't the problem either, but the problem is that our product gets approved for the market requires that each line of code (means also the functions invoked) is documented and tested. So the options are: build the math stuff from scratch and proof that the tests of them are justifing. or use the function of some one else and proof that its tests are justifing. So it woudl be more easy to write an own log as proofing some one else log works flawless, wouldn't it? $\endgroup$ – Zaibis Sep 16 '15 at 14:40
  • $\begingroup$ What kind of proofs are you looking for? Verifcation of correctness? Possible security flaws? Something else? The crlibm I liked to are verifyably correct. $\endgroup$ – mrf Sep 16 '15 at 14:41
  • $\begingroup$ so In other words, I'm looking for a simplification of a log while I can accept a certain inaccuracy $\endgroup$ – Zaibis Sep 16 '15 at 14:42
  • $\begingroup$ It's really hard to understand why your (most likely inferior) implementation would be preferable. To me it seems like testing (in the computer science sense) of a mathematical function is to make sure that the function is accurate (and not crashing on any input). Again, the library I linked satisfies those criteria. $\endgroup$ – mrf Sep 16 '15 at 14:45
  • $\begingroup$ I'm not completly aware of the kind of proof. But it is the advice of my boss that we don't use external sources because that would simply waste resources in the testing phase. so even if I would get convinced that doesn't change the fact that I have to write it. $\endgroup$ – Zaibis Sep 16 '15 at 14:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.