From 912ad2d3455228f3430bdf095536934ef42d1752 Mon Sep 17 00:00:00 2001 From: David Alessio Date: Thu, 14 Jul 2016 20:15:37 -0600 Subject: [PATCH] Fix bug in logf() algorithm that caused erroneous INFINITY results. --- libc/math/lib_erf.c | 3 ++- libc/math/lib_logf.c | 14 ++++++++------ 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/libc/math/lib_erf.c b/libc/math/lib_erf.c index 1923638dba..215c101720 100644 --- a/libc/math/lib_erf.c +++ b/libc/math/lib_erf.c @@ -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 diff --git a/libc/math/lib_logf.c b/libc/math/lib_logf.c index 2d097d228e..9abc5eff98 100644 --- a/libc/math/lib_logf.c +++ b/libc/math/lib_logf.c @@ -32,6 +32,8 @@ #include #include +#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; }