sergiotarxz/nuttx
sergiotarxz/nuttx
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

libc/math: fix log and logf calculations on ARMv7 (and maybe others)

Browse Source 
Probably this is a bug of a GCC, but on AMRv7 the code "if (relax_factor > 1)"
generates "bne.n" instruction if "relax_factor" is "int". Few lines above
"relax_factor *= LOG_RELAX_MULTIPLIER;" is done without overflow check hence
at some moment overflow occurs and "relax_factor" becomes a zero and condition
"if (relax_factor > 1)" becomes always evaluated to "true" hence "epsilon"
becomes zero always.

Probably this is not the best way to fix the bug (The best way is to report it
to GCC), but this change allows to get correct behavior of "log" and "logf" for
ARMv7 based MCUs

Signed-off-by: Petro Karashchenko <petro.karashchenko@gmail.com>
This commit is contained in:
Petro Karashchenko 2021-12-30 19:08:47 +02:00 committed by Xiang Xiao
parent 75aec04330
commit b859f2750b
3 changed files with 37 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
libs/libc/math/lib_log.c
View File

@ -39,13 +39,15 @@
* Pre-processor Definitions
****************************************************************************/
/* To avoid looping forever in particular corner cases, every LOGF_MAX_ITER
* the error criteria is relaxed by a factor LOGF_RELAX_MULTIPLIER.
#define DBL_MAX_EXP_X 700.0
/* To avoid looping forever in particular corner cases, every LOG_MAX_ITER
* the error criteria is relaxed by a factor LOG_RELAX_MULTIPLIER.
* todo: might need to adjust the double floating point version too.
*/
#define LOGF_MAX_ITER 10
#define LOGF_RELAX_MULTIPLIER 2
#define LOG_MAX_ITER 10
#define LOG_RELAX_MULTIPLIER 2.0
/****************************************************************************
* Public Functions
@ -62,7 +64,7 @@ double log(double x)
double y_old;
double ey;
double epsilon;
int relax_factor;
double relax_factor;
int iter;
y = 0.0;
@ -70,7 +72,7 @@ double log(double x)
epsilon = DBL_EPSILON;
iter = 0;
relax_factor = 1;
relax_factor = 1.0;
while (y > y_old + epsilon || y < y_old - epsilon)
{
@ -78,36 +80,36 @@ double log(double x)
ey = exp(y);
y -= (ey - x) / ey;
if (y > 700.0)
if (y > DBL_MAX_EXP_X)
{
y = 700.0;
y = DBL_MAX_EXP_X;
}
if (y < -700.0)
if (y < -DBL_MAX_EXP_X)
{
y = -700.0;
y = -DBL_MAX_EXP_X;
}
epsilon = (fabs(y) > 1.0) ? fabs(y) * DBL_EPSILON : DBL_EPSILON;
if (++iter >= LOGF_MAX_ITER)
if (++iter >= LOG_MAX_ITER)
{
relax_factor *= LOGF_RELAX_MULTIPLIER;
relax_factor *= LOG_RELAX_MULTIPLIER;
iter = 0;
}
if (relax_factor > 1)
if (relax_factor > 1.0)
{
epsilon *= relax_factor;
}
}
if (y == 700.0)
if (y == DBL_MAX_EXP_X)
{
return INFINITY;
}
if (y == -700.0)
if (y == -DBL_MAX_EXP_X)
{
return INFINITY;
}

8
libs/libc/math/lib_logf.c
View File

@ -44,7 +44,7 @@
*/
#define LOGF_MAX_ITER 10
#define LOGF_RELAX_MULTIPLIER 2
#define LOGF_RELAX_MULTIPLIER 2.0F
/****************************************************************************
* Public Functions
@ -60,7 +60,7 @@ float logf(float x)
float y_old;
float ey;
float epsilon;
int relax_factor;
float relax_factor;
int iter;
y = 0.0F;
@ -68,7 +68,7 @@ float logf(float x)
epsilon = FLT_EPSILON;
iter = 0;
relax_factor = 1;
relax_factor = 1.0F;
while (y > y_old + epsilon || y < y_old - epsilon)
{
@ -94,7 +94,7 @@ float logf(float x)
iter = 0;
}
if (relax_factor > 1)
if (relax_factor > 1.0F)
{
epsilon *= relax_factor;
}

24
libs/libc/math/lib_logl.c
View File

@ -35,6 +35,12 @@
#include <math.h>
#include <float.h>
/****************************************************************************
* Pre-processor Definitions
****************************************************************************/
#define LDBL_MAX_EXP_X 700.0
/****************************************************************************
* Public Functions
****************************************************************************/
@ -47,8 +53,8 @@ long double logl(long double x)
long double ey;
long double epsilon;
y = 0.0;
y_old = 1.0;
y = 0.0;
y_old = 1.0;
epsilon = LDBL_EPSILON;
while (y > y_old + epsilon || y < y_old - epsilon)
@ -57,23 +63,25 @@ long double logl(long double x)
ey = expl(y);
y -= (ey - x) / ey;
if (y > 700.0)
if (y > LDBL_MAX_EXP_X)
{
y = 700.0;
y = LDBL_MAX_EXP_X;
}
if (y < -700.0)
if (y < -LDBL_MAX_EXP_X)
{
y = -700.0;
y = -LDBL_MAX_EXP_X;
}
epsilon = (fabsl(y) > 1.0) ? fabsl(y) * LDBL_EPSILON : LDBL_EPSILON;
}
if (y == 700.0)
if (y == LDBL_MAX_EXP_X)
{
return INFINITY;
}
if (y == -700.0)
if (y == -LDBL_MAX_EXP_X)
{
return INFINITY;
}