my_log1p.c

#include <math.h>
#include <stdio.h>
#include <Python.h>
#include "numpy/npy_math.h"

/* Comparing numpy's fallback implementation of log1p with that of the C
   standard library.

   Compile with

   gcc -I$PYTHONBASE/lib/python2.7/site-packages/numpy/core/include \
       -I$PYTHONBASE/include/python2.7 \
       -L$PYTHONBASE/lib/python2.7/site-packages/numpy/core/lib \
       my_log1p.c -o my_log1p -lnpymath

   where PYTHONBASE points to a Python installation with numpy installed, e.g.
   PYTHONBASE=$(dirname $(dirname $(which python))) */


/* Fallback code path from numpy */
double my_npy_log1p(double x)
{
    if (npy_isinf(x) && x > 0) {
        return x;
    }
    else {
        const double u = 1. + x;
        const double d = u - 1.;

        if (d == 0) {
            return x;
        } else {
            return npy_log(u) * x / d;
        }
    }
}

int main()
{
    double x = 1.7976931348622732e+308;
    printf("Fallback: %f\n", my_npy_log1p(x));
    printf("The real thing: %f\n", npy_log1p(x));
    printf("libc's log1p: %f\n", log1p(x));
    return 0;
}