numeric stability of exp function

[lgarithm]

The logistic function yields slightly different results when defined in C and C++ source files:

static inline float logistic(float x) { return 1. / (1. + exp(-x)); }

[lgarithm]

#include <math.h>
#include <stdint.h>
#include <stdio.h>

static inline float logistic(float x) { return 1. / (1. + exp(-x)); }

void phex(float x)
{
    uint32_t *y = (uint32_t *)&x;
    printf("%08x", *y);
}

void test_logistic(float x)
{
    float y = logistic(x);
    printf("logistic(%f)=%f\n", x, y);
    phex(x);
    putchar(' ');
    phex(y);
    putchar('\n');
}

int main()
{
    float x;
    uint32_t *y = (uint32_t *)&x;
    *y = 0x3f788a4c;

    test_logistic(x);

    return 0;
}

compile as C and C++

gcc main.c   -o main_c
g++ main.cpp -o main_cpp

yields different results:

++ ./main_c
logistic(0.970860)=0.725291
3f788a4c 3f39acab
++ ./main_cpp
logistic(0.970860)=0.725291
3f788a4c 3f39acaa

[lgarithm]

softmax also have the same problem, so root cause is very likely to be exp.

[lgarithm]

In C, the argument type of exp is double: http://en.cppreference.com/w/c/numeric/math/exp