Fixes to asinh(), atanh(), and sinh(): The 'basic' expansions all exhibited bad cancellation errors near zero (<= 1E-10). This can be easily seen e.g. with x = 1E-30, the results are all zero though they should be extremely close to x. The cutoff values (1E-5, 1E-9) are chosen so that the next term in the Taylor series is negligible (for double). Functions could maybe be optimized to use only first term (x) and a smaller cutoff, just bigger than where the cancellation occurs.
This commit is contained in:
parent
018cbc75dc
commit
1025c3fb21
@ -49,6 +49,20 @@
|
||||
#ifdef CONFIG_HAVE_DOUBLE
|
||||
double asinh(double x)
|
||||
{
|
||||
return log(x + sqrt(x * x + 1));
|
||||
double y;
|
||||
double z = x * x;
|
||||
|
||||
if (z < 1E-9)
|
||||
{
|
||||
/* x - 1/6 * x^3 + 3/40 * x^5 - 5/112 * x^7 + ... */
|
||||
|
||||
y = x * (1 - z / 6.0);
|
||||
}
|
||||
else
|
||||
{
|
||||
y = log(x + sqrt(z + 1));
|
||||
}
|
||||
|
||||
return y;
|
||||
}
|
||||
#endif
|
||||
|
@ -49,6 +49,21 @@
|
||||
#ifdef CONFIG_HAVE_DOUBLE
|
||||
double atanh(double x)
|
||||
{
|
||||
return 0.5 * log((1 + x) / (1 - x));
|
||||
double y;
|
||||
|
||||
if (fabs(x) < 1E-5)
|
||||
{
|
||||
double z = x * x;
|
||||
|
||||
/* x + 1/3 * x^3 + 1/5 * x^5 + 1/7 * x^7 + ... */
|
||||
|
||||
y = x * (1 + z / 3.0);
|
||||
}
|
||||
else
|
||||
{
|
||||
y = log((1 + x) / (1 - x)) / 2.0;
|
||||
}
|
||||
|
||||
return y;
|
||||
}
|
||||
#endif
|
||||
|
@ -3,7 +3,7 @@
|
||||
*
|
||||
* This file is a part of NuttX:
|
||||
*
|
||||
* Copyright (C) 2012 Gregory Nutt. All rights reserved.
|
||||
* Copyright (C) 2012, 2015 Gregory Nutt. All rights reserved.
|
||||
* Ported by: Darcy Gong
|
||||
*
|
||||
* It derives from the Rhombs OS math library by Nick Johnson which has
|
||||
@ -41,7 +41,22 @@
|
||||
#ifdef CONFIG_HAVE_DOUBLE
|
||||
double sinh(double x)
|
||||
{
|
||||
x = exp(x);
|
||||
return ((x - (1.0 / x)) / 2.0);
|
||||
double y;
|
||||
double z;
|
||||
|
||||
if (fabs(x) < 1E-5)
|
||||
{
|
||||
/* x + 1/3! * x^3 + 1/5! * x^5 + 1/7! * x^7 + ... */
|
||||
|
||||
z = x * x;
|
||||
y = x * (1.0 + z / 6.0);
|
||||
}
|
||||
else
|
||||
{
|
||||
z = exp(x);
|
||||
y = (z - 1.0 / z) / 2.0;
|
||||
}
|
||||
|
||||
return y;
|
||||
}
|
||||
#endif
|
||||
|
Loading…
Reference in New Issue
Block a user