mirrors/phobos
1
0
Fork 0
mirror of https://github.com/dlang/phobos.git synced 2025-04-28 14:10:30 +03:00
Code Issues Projects Releases Packages Wiki Activity Actions

fix real precision use for inverse and hyperbolic sinh, tanh, scalbn, asinh, acosh (#7651)

Browse source
This commit is contained in:
Nicholas Wilson 2020-10-08 16:46:36 +08:00 • committed by GitHub
parent cf9b0c7b0d
commit f85ca8dbe8
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
1 changed files with 91 additions and 67 deletions
Download patch file Download diff file Expand all files Collapse all files

158
std/math.d
View file

@ -1322,39 +1322,45 @@ float cosh(float x) @safe pure nothrow @nogc { return cosh(cast(real) x); }
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
* )
*/
real sinh(real x) @safe pure nothrow @nogc
real sinh(real x) @safe pure nothrow @nogc { return _sinh(x); }
/// ditto
double sinh(double x) @safe pure nothrow @nogc { return _sinh(x); }
/// ditto
float sinh(float x) @safe pure nothrow @nogc { return _sinh(x); }
///
@safe unittest
{
enum sinh1 = (E - 1.0 / E) / 2;
import std.meta : AliasSeq;
static foreach (F; AliasSeq!(float, double, real))
{
assert(isIdentical(sinh(F(0.0)), F(0.0)));
assert(sinh(F(1.0)).approxEqual(F(sinh1)));
}
}
private F _sinh(F)(F x)
{
// sinh(x) = (exp(x)-exp(-x))/2;
// Very large arguments could cause an overflow, but
// the maximum value of x for which exp(x) + exp(-x)) != exp(x)
// is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80.
if (fabs(x) > real.mant_dig * LN2)
if (fabs(x) > F.mant_dig * F(LN2))
{
return copysign(0.5 * exp(fabs(x)), x);
return copysign(F(0.5) * exp(fabs(x)), x);
}
const real y = expm1(x);
return 0.5 * y / (y+1) * (y+2);
}
/// ditto
double sinh(double x) @safe pure nothrow @nogc { return sinh(cast(real) x); }
/// ditto
float sinh(float x) @safe pure nothrow @nogc { return sinh(cast(real) x); }
///
@safe unittest
{
assert(isIdentical(sinh(0.0), 0.0));
assert(sinh(1.0).approxEqual((E - 1.0 / E) / 2));
const y = expm1(x);
return F(0.5) * y / (y+1) * (y+2);
}
@safe @nogc nothrow unittest
{
assert(equalsDigit(sinh(1.0), (E - 1.0 / E) / 2, useDigits));
assert(equalsDigit(sinh(1.0L), real((E - 1.0 / E) / 2), useDigits));
}
/***********************************
* Calculates the hyperbolic tangent of x.
*
@ -1364,23 +1370,13 @@ float sinh(float x) @safe pure nothrow @nogc { return sinh(cast(real) x); }
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
* )
*/
real tanh(real x) @safe pure nothrow @nogc
{
// tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
if (fabs(x) > real.mant_dig * LN2)
{
return copysign(1, x);
}
const real y = expm1(2*x);
return y / (y + 2);
}
real tanh(real x) @safe pure nothrow @nogc { return _tanh(x); }
/// ditto
double tanh(double x) @safe pure nothrow @nogc { return tanh(cast(real) x); }
double tanh(double x) @safe pure nothrow @nogc { return _tanh(x); }
/// ditto
float tanh(float x) @safe pure nothrow @nogc { return tanh(cast(real) x); }
float tanh(float x) @safe pure nothrow @nogc { return _tanh(x); }
///
@safe unittest
@ -1389,9 +1385,21 @@ float tanh(float x) @safe pure nothrow @nogc { return tanh(cast(real) x); }
assert(tanh(1.0).approxEqual(sinh(1.0) / cosh(1.0)));
}
private F _tanh(F)(F x)
{
// tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
if (fabs(x) > F.mant_dig * F(LN2))
{
return copysign(1, x);
}
const y = expm1(2*x);
return y / (y + 2);
}
@safe @nogc nothrow unittest
{
assert(equalsDigit(tanh(1.0), sinh(1.0) / cosh(1.0), 15));
assert(equalsDigit(tanh(1.0L), sinh(1.0L) / cosh(1.0L), 15));
}
/***********************************
@ -1412,19 +1420,13 @@ float tanh(float x) @safe pure nothrow @nogc { return tanh(cast(real) x); }
* $(SV +$(INFIN),+$(INFIN))
* )
*/
real acosh(real x) @safe pure nothrow @nogc
{
if (x > 1/real.epsilon)
return LN2 + log(x);
else
return log(x + sqrt(x*x - 1));
}
real acosh(real x) @safe pure nothrow @nogc { return _acosh(x); }
/// ditto
double acosh(double x) @safe pure nothrow @nogc { return acosh(cast(real) x); }
double acosh(double x) @safe pure nothrow @nogc { return _acosh(x); }
/// ditto
float acosh(float x) @safe pure nothrow @nogc { return acosh(cast(real) x); }
float acosh(float x) @safe pure nothrow @nogc { return _acosh(x); }
///
@safe @nogc nothrow unittest
@ -1436,9 +1438,17 @@ float acosh(float x) @safe pure nothrow @nogc { return acosh(cast(real) x); }
assert(isNaN(acosh(0.5)));
}
private F _acosh(F)(F x) @safe pure nothrow @nogc
{
if (x > 1/F.epsilon)
return F(LN2) + log(x);
else
return log(x + sqrt(x*x - 1));
}
@safe @nogc nothrow unittest
{
assert(equalsDigit(acosh(cosh(3.0)), 3, useDigits));
assert(equalsDigit(acosh(cosh(3.0L)), 3.0L, useDigits));
}
/***********************************
@ -1457,20 +1467,13 @@ float acosh(float x) @safe pure nothrow @nogc { return acosh(cast(real) x); }
* $(SV $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
* )
*/
real asinh(real x) @safe pure nothrow @nogc
{
return (fabs(x) > 1 / real.epsilon)
// beyond this point, x*x + 1 == x*x
? copysign(LN2 + log(fabs(x)), x)
// sqrt(x*x + 1) == 1 + x * x / ( 1 + sqrt(x*x + 1) )
: copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
}
real asinh(real x) @safe pure nothrow @nogc { return _asinh(x); }
/// ditto
double asinh(double x) @safe pure nothrow @nogc { return asinh(cast(real) x); }
double asinh(double x) @safe pure nothrow @nogc { return _asinh(x); }
/// ditto
float asinh(float x) @safe pure nothrow @nogc { return asinh(cast(real) x); }
float asinh(float x) @safe pure nothrow @nogc { return _asinh(x); }
///
@safe @nogc nothrow unittest
@ -1482,9 +1485,18 @@ float asinh(float x) @safe pure nothrow @nogc { return asinh(cast(real) x); }
assert(isNaN(asinh(real.nan)));
}
private F _asinh(F)(F x)
{
return (fabs(x) > 1 / F.epsilon)
// beyond this point, x*x + 1 == x*x
? copysign(F(LN2) + log(fabs(x)), x)
// sqrt(x*x + 1) == 1 + x * x / ( 1 + sqrt(x*x + 1) )
: copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
}
@safe unittest
{
assert(equalsDigit(asinh(sinh(3.0)), 3, useDigits));
assert(equalsDigit(asinh(sinh(3.0L)), 3.0L, useDigits));
}
/***********************************
@ -4010,35 +4022,47 @@ real modf(real x, ref real i) @trusted nothrow @nogc
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
* )
*/
real scalbn(real x, int n) @safe pure nothrow @nogc
real scalbn(real x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
/// ditto
double scalbn(double x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
/// ditto
float scalbn(float x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
///
@safe pure nothrow @nogc unittest
{
assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
assert(scalbn(-real.infinity, 5) == -real.infinity);
assert(scalbn(2.0,10) == 2048.0);
assert(scalbn(2048.0f,-10) == 2.0f);
}
private F _scalbn(F)(F x, int n)
{
pragma(inline, true);
if (__ctfe)
{
// Handle special cases.
if (x == 0.0 || isInfinity(x))
if (x == F(0.0) || isInfinity(x))
return x;
}
return core.math.ldexp(x, n);
}
///
@safe pure nothrow @nogc unittest
{
assert(scalbn(0x1.2345678abcdefp0, 999) == 0x1.2345678abcdefp999);
assert(scalbn(-real.infinity, 5) == -real.infinity);
}
@safe pure nothrow @nogc unittest
{
// CTFE-able test
static assert(scalbn(0x1.2345678abcdefp0, 999) == 0x1.2345678abcdefp999);
static assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
static assert(scalbn(-real.infinity, 5) == -real.infinity);
// Test with large exponent delta n where the result is in bounds but 2.0L ^^ n is not.
enum initialExponent = real.min_exp + 2, resultExponent = real.max_exp - 2;
enum int n = resultExponent - initialExponent;
enum real x = 0x1.2345678abcdefp0 * (2.0L ^^ initialExponent);
enum real x = 0x1.2345678abcdefp0L * (2.0L ^^ initialExponent);
enum staticResult = scalbn(x, n);
static assert(staticResult == 0x1.2345678abcdefp0 * (2.0L ^^ resultExponent));
static assert(staticResult == 0x1.2345678abcdefp0L * (2.0L ^^ resultExponent));
assert(scalbn(x, n) == staticResult);
}