Fix bug in logf() algorithm that caused erroneous INFINITY results.

This commit is contained in:
David Alessio 2016-07-14 20:15:37 -06:00 committed by Gregory Nutt
parent 18059d6821
commit 912ad2d345
2 changed files with 10 additions and 7 deletions

View File

@ -67,6 +67,7 @@ double erf(double x)
sign = (x >= 0 ? 1 : -1);
t = 1.0/(1.0 + p*x);
return sign * (1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * expf(-x * x));
return sign * (1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t *
(double)expf(-x * x));
}
#endif

View File

@ -32,6 +32,8 @@
#include <math.h>
#include <float.h>
#define FLT_MAX_EXP_X 88.0F
/****************************************************************************
* Public Functions
****************************************************************************/
@ -50,25 +52,25 @@ float logf(float x)
ey = exp(y);
y -= (ey - x) / ey;
if (y > 700.0F)
if (y > FLT_MAX_EXP_X)
{
y = 700.0F;
y = FLT_MAX_EXP_X;
}
if (y < -700.0F)
if (y < -FLT_MAX_EXP_X)
{
y = -700.0F;
y = -FLT_MAX_EXP_X;
}
epsilon = (fabsf(y) > 1.0F) ? fabsf(y) * FLT_EPSILON : FLT_EPSILON;
}
if (y == 700.0F)
if (y == FLT_MAX_EXP_X)
{
return INFINITY;
}
if (y == -700.0F)
if (y == -FLT_MAX_EXP_X)
{
return INFINITY;
}