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

libs/libc/math: Small change to log() and logf() that improves accuracy and convergence time

Browse Source
This commit is contained in:
Rajan Gill 2018-09-16 15:14:15 -06:00 committed by Gregory Nutt
parent ed5620897f
commit d647127d14
2 changed files with 22 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -60,24 +60,24 @@ double log(double x)
{ {
double y; double y;
double y_old; double y_old;
double ey; double ney;
double epsilon; double epsilon;
int relax_factor; int relax_factor;
int i; int iter;
y = 0.0; y = 0.0;
y_old = 1.0; y_old = 1.0;
epsilon = DBL_EPSILON; epsilon = DBL_EPSILON;
i = 0; iter = 0;
relax_factor = 1; relax_factor = 1;
while (y > y_old + epsilon || y < y_old - epsilon) while (y > y_old + epsilon || y < y_old - epsilon)
{ {
y_old = y; y_old = y;
ey = exp(y); ney = exp(-y);
y -= (ey - x) / ey; y -= 1 - x * ney;
if (y > 700.0) if (y > 700.0)
{ {
@ -91,13 +91,16 @@ double log(double x)
epsilon = (fabs(y) > 1.0) ? fabs(y) * DBL_EPSILON : DBL_EPSILON; epsilon = (fabs(y) > 1.0) ? fabs(y) * DBL_EPSILON : DBL_EPSILON;
if (++i >= LOGF_MAX_ITER) if (++iter >= LOGF_MAX_ITER)
{ {
relax_factor *= LOGF_RELAX_MULTIPLIER; relax_factor *= LOGF_RELAX_MULTIPLIER;
i = 0; iter = 0;
} }
epsilon *= relax_factor; if (relax_factor > 1)
{
epsilon *= relax_factor;
}
} }
if (y == 700.0) if (y == 700.0)

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

@ -58,23 +58,23 @@ float logf(float x)
{ {
float y; float y;
float y_old; float y_old;
float ey; float ney;
float epsilon; float epsilon;
int relax_factor; int relax_factor;
int i; int iter;
y = 0.0F; y = 0.0F;
y_old = 1.0F; y_old = 1.0F;
epsilon = FLT_EPSILON; epsilon = FLT_EPSILON;
i = 0; iter = 0;
relax_factor = 1; relax_factor = 1;
while (y > y_old + epsilon || y < y_old - epsilon) while (y > y_old + epsilon || y < y_old - epsilon)
{ {
y_old = y; y_old = y;
ey = expf(y); ney = expf(-y);
y -= (ey - x) / ey; y -= 1 - x * ney;
if (y > FLT_MAX_EXP_X) if (y > FLT_MAX_EXP_X)
{ {
@ -88,13 +88,16 @@ float logf(float x)
epsilon = (fabsf(y) > 1.0F) ? fabsf(y) * FLT_EPSILON : FLT_EPSILON; epsilon = (fabsf(y) > 1.0F) ? fabsf(y) * FLT_EPSILON : FLT_EPSILON;
if (++i >= LOGF_MAX_ITER) if (++iter >= LOGF_MAX_ITER)
{ {
relax_factor *= LOGF_RELAX_MULTIPLIER; relax_factor *= LOGF_RELAX_MULTIPLIER;
i = 0; iter = 0;
} }
epsilon *= relax_factor; if (relax_factor > 1)
{
epsilon *= relax_factor;
}
} }
if (y == FLT_MAX_EXP_X) if (y == FLT_MAX_EXP_X)