Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
6833 lines (6234 sloc) 203.563 kB
// Written in the D programming language.
/**
* Contains the elementary mathematical functions (powers, roots,
* and trigonometric functions), and low-level floating-point operations.
* Mathematical special functions are available in $(D std.mathspecial).
*
$(SCRIPT inhibitQuickIndex = 1;)
$(DIVC quickindex,
$(BOOKTABLE ,
$(TR $(TH Category) $(TH Members) )
$(TR $(TDNW Constants) $(TD
$(MYREF E) $(MYREF PI) $(MYREF PI_2) $(MYREF PI_4) $(MYREF M_1_PI)
$(MYREF M_2_PI) $(MYREF M_2_SQRTPI) $(MYREF LN10) $(MYREF LN2)
$(MYREF LOG2) $(MYREF LOG2E) $(MYREF LOG2T) $(MYREF LOG10E)
$(MYREF SQRT2) $(MYREF SQRT1_2)
))
$(TR $(TDNW Classics) $(TD
$(MYREF abs) $(MYREF fabs) $(MYREF sqrt) $(MYREF cbrt) $(MYREF hypot) $(MYREF poly)
))
$(TR $(TDNW Trigonometry) $(TD
$(MYREF sin) $(MYREF cos) $(MYREF tan) $(MYREF asin) $(MYREF acos)
$(MYREF atan) $(MYREF atan2) $(MYREF sinh) $(MYREF cosh) $(MYREF tanh)
$(MYREF asinh) $(MYREF acosh) $(MYREF atanh) $(MYREF expi)
))
$(TR $(TDNW Rounding) $(TD
$(MYREF ceil) $(MYREF floor) $(MYREF round) $(MYREF lround)
$(MYREF trunc) $(MYREF rint) $(MYREF lrint) $(MYREF nearbyint)
$(MYREF rndtol)
))
$(TR $(TDNW Exponentiation & Logarithms) $(TD
$(MYREF pow) $(MYREF exp) $(MYREF exp2) $(MYREF expm1) $(MYREF ldexp)
$(MYREF frexp) $(MYREF log) $(MYREF log2) $(MYREF log10) $(MYREF logb)
$(MYREF ilogb) $(MYREF log1p) $(MYREF scalbn)
))
$(TR $(TDNW Modulus) $(TD
$(MYREF fmod) $(MYREF modf) $(MYREF remainder)
))
$(TR $(TDNW Floating-point operations) $(TD
$(MYREF approxEqual) $(MYREF feqrel) $(MYREF fdim) $(MYREF fmax)
$(MYREF fmin) $(MYREF fma) $(MYREF nextDown) $(MYREF nextUp)
$(MYREF nextafter) $(MYREF NaN) $(MYREF getNaNPayload)
))
$(TR $(TDNW Introspection) $(TD
$(MYREF isFinite) $(MYREF isIdentical) $(MYREF isInfinity) $(MYREF isNaN)
$(MYREF isNormal) $(MYREF isSubnormal) $(MYREF signbit) $(MYREF sgn)
$(MYREF copysign)
))
$(TR $(TDNW Complex Numbers) $(TD
$(MYREF abs) $(MYREF conj) $(MYREF sin) $(MYREF cos) $(MYREF expi)
))
$(TR $(TDNW Hardware Control) $(TD
$(MYREF IeeeFlags) $(MYREF FloatingPointControl)
))
)
)
* The functionality closely follows the IEEE754-2008 standard for
* floating-point arithmetic, including the use of camelCase names rather
* than C99-style lower case names. All of these functions behave correctly
* when presented with an infinity or NaN.
*
* The following IEEE 'real' formats are currently supported:
* $(UL
* $(LI 64 bit Big-endian 'double' (eg PowerPC))
* $(LI 128 bit Big-endian 'quadruple' (eg SPARC))
* $(LI 64 bit Little-endian 'double' (eg x86-SSE2))
* $(LI 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium))
* $(LI 128 bit Little-endian 'quadruple' (not implemented on any known processor!))
* $(LI Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support)
* )
* Unlike C, there is no global 'errno' variable. Consequently, almost all of
* these functions are pure nothrow.
*
* Status:
* The semantics and names of feqrel and approxEqual will be revised.
*
* Macros:
* WIKI = Phobos/StdMath
*
* TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
* <caption>Special Values</caption>
* $0</table>
* SVH = $(TR $(TH $1) $(TH $2))
* SV = $(TR $(TD $1) $(TD $2))
* TH3 = $(TR $(TH $1) $(TH $2) $(TH $3))
* TD3 = $(TR $(TD $1) $(TD $2) $(TD $3))
*
* NAN = $(RED NAN)
* SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
* GAMMA = &#915;
* THETA = &theta;
* INTEGRAL = &#8747;
* INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
* POWER = $1<sup>$2</sup>
* SUB = $1<sub>$2</sub>
* BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
* CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
* PLUSMN = &plusmn;
* INFIN = &infin;
* PLUSMNINF = &plusmn;&infin;
* PI = &pi;
* LT = &lt;
* GT = &gt;
* SQRT = &radic;
* HALF = &frac12;
*
* Copyright: Copyright Digital Mars 2000 - 2011.
* D implementations of tan, atan, atan2, exp, expm1, exp2, log, log10, log1p,
* log2, floor, ceil and lrint functions are based on the CEPHES math library,
* which is Copyright (C) 2001 Stephen L. Moshier <steve@moshier.net>
* and are incorporated herein by permission of the author. The author
* reserves the right to distribute this material elsewhere under different
* copying permissions. These modifications are distributed here under
* the following terms:
* License: $(WEB www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
* Authors: $(WEB digitalmars.com, Walter Bright), Don Clugston,
* Conversion of CEPHES math library to D by Iain Buclaw
* Source: $(PHOBOSSRC std/_math.d)
*/
module std.math;
version (Win64)
{
version (D_InlineAsm_X86_64)
version = Win64_DMD_InlineAsm;
}
import core.stdc.math;
import std.traits;
version(LDC)
{
import ldc.intrinsics;
}
version(DigitalMars)
{
version = INLINE_YL2X; // x87 has opcodes for these
}
version (X86) version = X86_Any;
version (X86_64) version = X86_Any;
version (PPC) version = PPC_Any;
version (PPC64) version = PPC_Any;
version(D_InlineAsm_X86)
{
version = InlineAsm_X86_Any;
}
else version(D_InlineAsm_X86_64)
{
version = InlineAsm_X86_Any;
}
version(unittest)
{
import core.stdc.stdio;
static if(real.sizeof > double.sizeof)
enum uint useDigits = 16;
else
enum uint useDigits = 15;
/******************************************
* Compare floating point numbers to n decimal digits of precision.
* Returns:
* 1 match
* 0 nomatch
*/
private bool equalsDigit(real x, real y, uint ndigits)
{
if (signbit(x) != signbit(y))
return 0;
if (isInfinity(x) && isInfinity(y))
return 1;
if (isInfinity(x) || isInfinity(y))
return 0;
if (isNaN(x) && isNaN(y))
return 1;
if (isNaN(x) || isNaN(y))
return 0;
char[30] bufx;
char[30] bufy;
assert(ndigits < bufx.length);
int ix;
int iy;
version(CRuntime_Microsoft)
alias real_t = double;
else
alias real_t = real;
ix = sprintf(bufx.ptr, "%.*Lg", ndigits, cast(real_t) x);
iy = sprintf(bufy.ptr, "%.*Lg", ndigits, cast(real_t) y);
assert(ix < bufx.length && ix > 0);
assert(ix < bufy.length && ix > 0);
return bufx[0 .. ix] == bufy[0 .. iy];
}
}
private:
// The following IEEE 'real' formats are currently supported.
version(LittleEndian)
{
static assert(real.mant_dig == 53 || real.mant_dig == 64
|| real.mant_dig == 113,
"Only 64-bit, 80-bit, and 128-bit reals"~
" are supported for LittleEndian CPUs");
}
else
{
static assert(real.mant_dig == 53 || real.mant_dig == 106
|| real.mant_dig == 113,
"Only 64-bit and 128-bit reals are supported for BigEndian CPUs."~
" double-double reals have partial support");
}
// Underlying format exposed through floatTraits
enum RealFormat
{
ieeeHalf,
ieeeSingle,
ieeeDouble,
ieeeExtended, // x87 80-bit real
ieeeExtended53, // x87 real rounded to precision of double.
ibmExtended, // IBM 128-bit extended
ieeeQuadruple,
}
// Constants used for extracting the components of the representation.
// They supplement the built-in floating point properties.
template floatTraits(T)
{
// EXPMASK is a ushort mask to select the exponent portion (without sign)
// EXPSHIFT is the number of bits the exponent is left-shifted by in its ushort
// EXPPOS_SHORT is the index of the exponent when represented as a ushort array.
// SIGNPOS_BYTE is the index of the sign when represented as a ubyte array.
// RECIP_EPSILON is the value such that (smallest_subnormal) * RECIP_EPSILON == T.min_normal
enum T RECIP_EPSILON = (1/T.epsilon);
static if (T.mant_dig == 24)
{
// Single precision float
enum ushort EXPMASK = 0x7F80;
enum ushort EXPSHIFT = 7;
enum ushort EXPBIAS = 0x3F00;
enum uint EXPMASK_INT = 0x7F80_0000;
enum uint MANTISSAMASK_INT = 0x007F_FFFF;
enum realFormat = RealFormat.ieeeSingle;
version(LittleEndian)
{
enum EXPPOS_SHORT = 1;
enum SIGNPOS_BYTE = 3;
}
else
{
enum EXPPOS_SHORT = 0;
enum SIGNPOS_BYTE = 0;
}
}
else static if (T.mant_dig == 53)
{
static if (T.sizeof == 8)
{
// Double precision float, or real == double
enum ushort EXPMASK = 0x7FF0;
enum ushort EXPSHIFT = 4;
enum ushort EXPBIAS = 0x3FE0;
enum uint EXPMASK_INT = 0x7FF0_0000;
enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only
enum realFormat = RealFormat.ieeeDouble;
version(LittleEndian)
{
enum EXPPOS_SHORT = 3;
enum SIGNPOS_BYTE = 7;
}
else
{
enum EXPPOS_SHORT = 0;
enum SIGNPOS_BYTE = 0;
}
}
else static if (T.sizeof == 12)
{
// Intel extended real80 rounded to double
enum ushort EXPMASK = 0x7FFF;
enum ushort EXPSHIFT = 0;
enum ushort EXPBIAS = 0x3FFE;
enum realFormat = RealFormat.ieeeExtended53;
version(LittleEndian)
{
enum EXPPOS_SHORT = 4;
enum SIGNPOS_BYTE = 9;
}
else
{
enum EXPPOS_SHORT = 0;
enum SIGNPOS_BYTE = 0;
}
}
else
static assert(false, "No traits support for " ~ T.stringof);
}
else static if (T.mant_dig == 64)
{
// Intel extended real80
enum ushort EXPMASK = 0x7FFF;
enum ushort EXPSHIFT = 0;
enum ushort EXPBIAS = 0x3FFE;
enum realFormat = RealFormat.ieeeExtended;
version(LittleEndian)
{
enum EXPPOS_SHORT = 4;
enum SIGNPOS_BYTE = 9;
}
else
{
enum EXPPOS_SHORT = 0;
enum SIGNPOS_BYTE = 0;
}
}
else static if (T.mant_dig == 113)
{
// Quadruple precision float
enum ushort EXPMASK = 0x7FFF;
enum ushort EXPSHIFT = 0;
enum ushort EXPBIAS = 0x3FFF;
enum realFormat = RealFormat.ieeeQuadruple;
version(LittleEndian)
{
enum EXPPOS_SHORT = 7;
enum SIGNPOS_BYTE = 15;
}
else
{
enum EXPPOS_SHORT = 0;
enum SIGNPOS_BYTE = 0;
}
}
else static if (T.mant_dig == 106)
{
// IBM Extended doubledouble
enum ushort EXPMASK = 0x7FF0;
enum ushort EXPSHIFT = 4;
enum realFormat = RealFormat.ibmExtended;
// the exponent byte is not unique
version(LittleEndian)
{
enum EXPPOS_SHORT = 7; // [3] is also an exp short
enum SIGNPOS_BYTE = 15;
}
else
{
enum EXPPOS_SHORT = 0; // [4] is also an exp short
enum SIGNPOS_BYTE = 0;
}
}
else
static assert(false, "No traits support for " ~ T.stringof);
}
// These apply to all floating-point types
version(LittleEndian)
{
enum MANTISSA_LSB = 0;
enum MANTISSA_MSB = 1;
}
else
{
enum MANTISSA_LSB = 1;
enum MANTISSA_MSB = 0;
}
// Common code for math implementations.
// Helper for floor/ceil
T floorImpl(T)(const T x) @trusted pure nothrow @nogc
{
alias F = floatTraits!(T);
// Take care not to trigger library calls from the compiler,
// while ensuring that we don't get defeated by some optimizers.
union floatBits
{
T rv;
ushort[T.sizeof/2] vu;
}
floatBits y = void;
y.rv = x;
// Find the exponent (power of 2)
static if (F.realFormat == RealFormat.ieeeSingle)
{
int exp = ((y.vu[F.EXPPOS_SHORT] >> 7) & 0xff) - 0x7f;
version (LittleEndian)
int pos = 0;
else
int pos = 3;
}
else static if (F.realFormat == RealFormat.ieeeDouble)
{
int exp = ((y.vu[F.EXPPOS_SHORT] >> 4) & 0x7ff) - 0x3ff;
version (LittleEndian)
int pos = 0;
else
int pos = 3;
}
else static if (F.realFormat == RealFormat.ieeeExtended)
{
int exp = (y.vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff;
version (LittleEndian)
int pos = 0;
else
int pos = 4;
}
else static if (F.realFormat == RealFormat.ieeeQuadruple)
{
int exp = (y.vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff;
version (LittleEndian)
int pos = 0;
else
int pos = 7;
}
else
static assert(false, "Not implemented for this architecture");
if (exp < 0)
{
if (x < 0.0)
return -1.0;
else
return 0.0;
}
exp = (T.mant_dig - 1) - exp;
// Zero 16 bits at a time.
while (exp >= 16)
{
version (LittleEndian)
y.vu[pos++] = 0;
else
y.vu[pos--] = 0;
exp -= 16;
}
// Clear the remaining bits.
if (exp > 0)
y.vu[pos] &= 0xffff ^ ((1 << exp) - 1);
if ((x < 0.0) && (x != y.rv))
y.rv -= 1.0;
return y.rv;
}
public:
// Values obtained from Wolfram Alpha. 116 bits ought to be enough for anybody.
// Wolfram Alpha LLC. 2011. Wolfram|Alpha. http://www.wolframalpha.com/input/?i=e+in+base+16 (access July 6, 2011).
enum real E = 0x1.5bf0a8b1457695355fb8ac404e7a8p+1L; /** e = 2.718281... */
enum real LOG2T = 0x1.a934f0979a3715fc9257edfe9b5fbp+1L; /** $(SUB log, 2)10 = 3.321928... */
enum real LOG2E = 0x1.71547652b82fe1777d0ffda0d23a8p+0L; /** $(SUB log, 2)e = 1.442695... */
enum real LOG2 = 0x1.34413509f79fef311f12b35816f92p-2L; /** $(SUB log, 10)2 = 0.301029... */
enum real LOG10E = 0x1.bcb7b1526e50e32a6ab7555f5a67cp-2L; /** $(SUB log, 10)e = 0.434294... */
enum real LN2 = 0x1.62e42fefa39ef35793c7673007e5fp-1L; /** ln 2 = 0.693147... */
enum real LN10 = 0x1.26bb1bbb5551582dd4adac5705a61p+1L; /** ln 10 = 2.302585... */
enum real PI = 0x1.921fb54442d18469898cc51701b84p+1L; /** $(_PI) = 3.141592... */
enum real PI_2 = PI/2; /** $(PI) / 2 = 1.570796... */
enum real PI_4 = PI/4; /** $(PI) / 4 = 0.785398... */
enum real M_1_PI = 0x1.45f306dc9c882a53f84eafa3ea69cp-2L; /** 1 / $(PI) = 0.318309... */
enum real M_2_PI = 2*M_1_PI; /** 2 / $(PI) = 0.636619... */
enum real M_2_SQRTPI = 0x1.20dd750429b6d11ae3a914fed7fd8p+0L; /** 2 / $(SQRT)$(PI) = 1.128379... */
enum real SQRT2 = 0x1.6a09e667f3bcc908b2fb1366ea958p+0L; /** $(SQRT)2 = 1.414213... */
enum real SQRT1_2 = SQRT2/2; /** $(SQRT)$(HALF) = 0.707106... */
// Note: Make sure the magic numbers in compiler backend for x87 match these.
/***********************************
* Calculates the absolute value of a number
*
* Params:
* Num = (template parameter) type of number
* x = real number value
* z = complex number value
* y = imaginary number value
*
* Returns:
* The absolute value of the number. If floating-point or integral,
* the return type will be the same as the input; if complex or
* imaginary, the returned value will be the corresponding floating
* point type.
*
* For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) )
* = hypot(z.re, z.im).
*/
Num abs(Num)(Num x) @safe pure nothrow
if (is(typeof(Num.init >= 0)) && is(typeof(-Num.init)) &&
!(is(Num* : const(ifloat*)) || is(Num* : const(idouble*))
|| is(Num* : const(ireal*))))
{
static if (isFloatingPoint!(Num))
return fabs(x);
else
return x>=0 ? x : -x;
}
/// ditto
auto abs(Num)(Num z) @safe pure nothrow @nogc
if (is(Num* : const(cfloat*)) || is(Num* : const(cdouble*))
|| is(Num* : const(creal*)))
{
return hypot(z.re, z.im);
}
/// ditto
auto abs(Num)(Num y) @safe pure nothrow @nogc
if (is(Num* : const(ifloat*)) || is(Num* : const(idouble*))
|| is(Num* : const(ireal*)))
{
return fabs(y.im);
}
/// ditto
@safe pure nothrow @nogc unittest
{
assert(isIdentical(abs(-0.0L), 0.0L));
assert(isNaN(abs(real.nan)));
assert(abs(-real.infinity) == real.infinity);
assert(abs(-3.2Li) == 3.2L);
assert(abs(71.6Li) == 71.6L);
assert(abs(-56) == 56);
assert(abs(2321312L) == 2321312L);
assert(abs(-1L+1i) == sqrt(2.0L));
}
@safe pure nothrow @nogc unittest
{
import std.typetuple;
foreach (T; TypeTuple!(float, double, real))
{
T f = 3;
assert(abs(f) == f);
assert(abs(-f) == f);
}
foreach (T; TypeTuple!(cfloat, cdouble, creal))
{
T f = -12+3i;
assert(abs(f) == hypot(f.re, f.im));
assert(abs(-f) == hypot(f.re, f.im));
}
}
/***********************************
* Complex conjugate
*
* conj(x + iy) = x - iy
*
* Note that z * conj(z) = $(POWER z.re, 2) - $(POWER z.im, 2)
* is always a real number
*/
auto conj(Num)(Num z) @safe pure nothrow @nogc
if (is(Num* : const(cfloat*)) || is(Num* : const(cdouble*))
|| is(Num* : const(creal*)))
{
//FIXME
//Issue 14206
static if(is(Num* : const(cdouble*)))
return cast(cdouble) conj(cast(creal)z);
else
return z.re - z.im*1fi;
}
/** ditto */
auto conj(Num)(Num y) @safe pure nothrow @nogc
if (is(Num* : const(ifloat*)) || is(Num* : const(idouble*))
|| is(Num* : const(ireal*)))
{
return -y;
}
///
@safe pure nothrow @nogc unittest
{
creal c = 7 + 3Li;
assert(conj(c) == 7-3Li);
ireal z = -3.2Li;
assert(conj(z) == -z);
}
//Issue 14206
@safe pure nothrow @nogc unittest
{
cdouble c = 7 + 3i;
assert(conj(c) == 7-3i);
idouble z = -3.2i;
assert(conj(z) == -z);
}
//Issue 14206
@safe pure nothrow @nogc unittest
{
cfloat c = 7f + 3fi;
assert(conj(c) == 7f-3fi);
ifloat z = -3.2fi;
assert(conj(z) == -z);
}
/***********************************
* Returns cosine of x. x is in radians.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH cos(x)) $(TH invalid?))
* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) )
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes) )
* )
* Bugs:
* Results are undefined if |x| >= $(POWER 2,64).
*/
real cos(real x) @safe pure nothrow @nogc; /* intrinsic */
//FIXME
///ditto
double cos(double x) @safe pure nothrow @nogc { return cos(cast(real)x); }
//FIXME
///ditto
float cos(float x) @safe pure nothrow @nogc { return cos(cast(real)x); }
/***********************************
* Returns $(WEB en.wikipedia.org/wiki/Sine, sine) of x. x is in $(WEB en.wikipedia.org/wiki/Radian, radians).
*
* $(TABLE_SV
* $(TH3 x , sin(x) , invalid?)
* $(TD3 $(NAN) , $(NAN) , yes )
* $(TD3 $(PLUSMN)0.0, $(PLUSMN)0.0, no )
* $(TD3 $(PLUSMNINF), $(NAN) , yes )
* )
*
* Params:
* x = angle in radians (not degrees)
* Returns:
* sine of x
* See_Also:
* $(MYREF cos), $(MYREF tan), $(MYREF asin)
* Bugs:
* Results are undefined if |x| >= $(POWER 2,64).
*/
real sin(real x) @safe pure nothrow @nogc; /* intrinsic */
//FIXME
///ditto
double sin(double x) @safe pure nothrow @nogc { return sin(cast(real)x); }
//FIXME
///ditto
float sin(float x) @safe pure nothrow @nogc { return sin(cast(real)x); }
///
unittest
{
import std.math : sin, PI;
import std.stdio : writefln;
void someFunc()
{
real x = 30.0;
auto result = sin(x * (PI / 180)); // convert degrees to radians
writefln("The sine of %s degrees is %s", x, result);
}
}
/***********************************
* Returns sine for complex and imaginary arguments.
*
* sin(z) = sin(z.re)*cosh(z.im) + cos(z.re)*sinh(z.im)i
*
* If both sin($(THETA)) and cos($(THETA)) are required,
* it is most efficient to use expi($(THETA)).
*/
creal sin(creal z) @safe pure nothrow @nogc
{
creal cs = expi(z.re);
creal csh = coshisinh(z.im);
return cs.im * csh.re + cs.re * csh.im * 1i;
}
/** ditto */
ireal sin(ireal y) @safe pure nothrow @nogc
{
return cosh(y.im)*1i;
}
///
@safe pure nothrow @nogc unittest
{
assert(sin(0.0+0.0i) == 0.0);
assert(sin(2.0+0.0i) == sin(2.0L) );
}
/***********************************
* cosine, complex and imaginary
*
* cos(z) = cos(z.re)*cosh(z.im) - sin(z.re)*sinh(z.im)i
*/
creal cos(creal z) @safe pure nothrow @nogc
{
creal cs = expi(z.re);
creal csh = coshisinh(z.im);
return cs.re * csh.re - cs.im * csh.im * 1i;
}
/** ditto */
real cos(ireal y) @safe pure nothrow @nogc
{
return cosh(y.im);
}
///
@safe pure nothrow @nogc unittest
{
assert(cos(0.0+0.0i)==1.0);
assert(cos(1.3L+0.0i)==cos(1.3L));
assert(cos(5.2Li)== cosh(5.2L));
}
/****************************************************************************
* Returns tangent of x. x is in radians.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH tan(x)) $(TH invalid?))
* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes))
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
* $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD yes))
* )
*/
real tan(real x) @trusted pure nothrow @nogc
{
version(D_InlineAsm_X86)
{
asm pure nothrow @nogc
{
fld x[EBP] ; // load theta
fxam ; // test for oddball values
fstsw AX ;
sahf ;
jc trigerr ; // x is NAN, infinity, or empty
// 387's can handle subnormals
SC18: fptan ;
fstp ST(0) ; // dump X, which is always 1
fstsw AX ;
sahf ;
jnp Lret ; // C2 = 1 (x is out of range)
// Do argument reduction to bring x into range
fldpi ;
fxch ;
SC17: fprem1 ;
fstsw AX ;
sahf ;
jp SC17 ;
fstp ST(1) ; // remove pi from stack
jmp SC18 ;
trigerr:
jnp Lret ; // if theta is NAN, return theta
fstp ST(0) ; // dump theta
}
return real.nan;
Lret: {}
}
else version(D_InlineAsm_X86_64)
{
version (Win64)
{
asm pure nothrow @nogc
{
fld real ptr [RCX] ; // load theta
}
}
else
{
asm pure nothrow @nogc
{
fld x[RBP] ; // load theta
}
}
asm pure nothrow @nogc
{
fxam ; // test for oddball values
fstsw AX ;
test AH,1 ;
jnz trigerr ; // x is NAN, infinity, or empty
// 387's can handle subnormals
SC18: fptan ;
fstp ST(0) ; // dump X, which is always 1
fstsw AX ;
test AH,4 ;
jz Lret ; // C2 = 1 (x is out of range)
// Do argument reduction to bring x into range
fldpi ;
fxch ;
SC17: fprem1 ;
fstsw AX ;
test AH,4 ;
jnz SC17 ;
fstp ST(1) ; // remove pi from stack
jmp SC18 ;
trigerr:
test AH,4 ;
jz Lret ; // if theta is NAN, return theta
fstp ST(0) ; // dump theta
}
return real.nan;
Lret: {}
}
else
{
// Coefficients for tan(x)
static immutable real[3] P = [
-1.7956525197648487798769E7L,
1.1535166483858741613983E6L,
-1.3093693918138377764608E4L,
];
static immutable real[5] Q = [
-5.3869575592945462988123E7L,
2.5008380182335791583922E7L,
-1.3208923444021096744731E6L,
1.3681296347069295467845E4L,
1.0000000000000000000000E0L,
];
// PI/4 split into three parts.
enum real P1 = 7.853981554508209228515625E-1L;
enum real P2 = 7.946627356147928367136046290398E-9L;
enum real P3 = 3.061616997868382943065164830688E-17L;
// Special cases.
if (x == 0.0 || isNaN(x))
return x;
if (isInfinity(x))
return real.nan;
// Make argument positive but save the sign.
bool sign = false;
if (signbit(x))
{
sign = true;
x = -x;
}
// Compute x mod PI/4.
real y = floor(x / PI_4);
// Strip high bits of integer part.
real z = ldexp(y, -4);
// Compute y - 16 * (y / 16).
z = y - ldexp(floor(z), 4);
// Integer and fraction part modulo one octant.
int j = cast(int)(z);
// Map zeros and singularities to origin.
if (j & 1)
{
j += 1;
y += 1.0;
}
z = ((x - y * P1) - y * P2) - y * P3;
real zz = z * z;
if (zz > 1.0e-20L)
y = z + z * (zz * poly(zz, P) / poly(zz, Q));
else
y = z;
if (j & 2)
y = -1.0 / y;
return (sign) ? -y : y;
}
}
@safe nothrow @nogc unittest
{
static real[2][] vals = // angle,tan
[
[ 0, 0],
[ .5, .5463024898],
[ 1, 1.557407725],
[ 1.5, 14.10141995],
[ 2, -2.185039863],
[ 2.5,-.7470222972],
[ 3, -.1425465431],
[ 3.5, .3745856402],
[ 4, 1.157821282],
[ 4.5, 4.637332055],
[ 5, -3.380515006],
[ 5.5,-.9955840522],
[ 6, -.2910061914],
[ 6.5, .2202772003],
[ 10, .6483608275],
// special angles
[ PI_4, 1],
//[ PI_2, real.infinity], // PI_2 is not _exactly_ pi/2.
[ 3*PI_4, -1],
[ PI, 0],
[ 5*PI_4, 1],
//[ 3*PI_2, -real.infinity],
[ 7*PI_4, -1],
[ 2*PI, 0],
];
int i;
for (i = 0; i < vals.length; i++)
{
real x = vals[i][0];
real r = vals[i][1];
real t = tan(x);
//printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
if (!isIdentical(r, t)) assert(fabs(r-t) <= .0000001);
x = -x;
r = -r;
t = tan(x);
//printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
if (!isIdentical(r, t) && !(r!=r && t!=t)) assert(fabs(r-t) <= .0000001);
}
// overflow
assert(isNaN(tan(real.infinity)));
assert(isNaN(tan(-real.infinity)));
// NaN propagation
assert(isIdentical( tan(NaN(0x0123L)), NaN(0x0123L) ));
}
unittest
{
assert(equalsDigit(tan(PI / 3), std.math.sqrt(3.0), useDigits));
}
/***************
* Calculates the arc cosine of x,
* returning a value ranging from 0 to $(PI).
*
* $(TABLE_SV
* $(TR $(TH x) $(TH acos(x)) $(TH invalid?))
* $(TR $(TD $(GT)1.0) $(TD $(NAN)) $(TD yes))
* $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD yes))
* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes))
* )
*/
real acos(real x) @safe pure nothrow @nogc
{
return atan2(sqrt(1-x*x), x);
}
/// ditto
double acos(double x) @safe pure nothrow @nogc { return acos(cast(real)x); }
/// ditto
float acos(float x) @safe pure nothrow @nogc { return acos(cast(real)x); }
unittest
{
assert(equalsDigit(acos(0.5), std.math.PI / 3, useDigits));
}
/***************
* Calculates the arc sine of x,
* returning a value ranging from -$(PI)/2 to $(PI)/2.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH asin(x)) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
* $(TR $(TD $(GT)1.0) $(TD $(NAN)) $(TD yes))
* $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD yes))
* )
*/
real asin(real x) @safe pure nothrow @nogc
{
return atan2(x, sqrt(1-x*x));
}
/// ditto
double asin(double x) @safe pure nothrow @nogc { return asin(cast(real)x); }
/// ditto
float asin(float x) @safe pure nothrow @nogc { return asin(cast(real)x); }
unittest
{
assert(equalsDigit(asin(0.5), PI / 6, useDigits));
}
/***************
* Calculates the arc tangent of x,
* returning a value ranging from -$(PI)/2 to $(PI)/2.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH atan(x)) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes))
* )
*/
real atan(real x) @safe pure nothrow @nogc
{
version(InlineAsm_X86_Any)
{
return atan2(x, 1.0L);
}
else
{
// Coefficients for atan(x)
static immutable real[5] P = [
-5.0894116899623603312185E1L,
-9.9988763777265819915721E1L,
-6.3976888655834347413154E1L,
-1.4683508633175792446076E1L,
-8.6863818178092187535440E-1L,
];
static immutable real[6] Q = [
1.5268235069887081006606E2L,
3.9157570175111990631099E2L,
3.6144079386152023162701E2L,
1.4399096122250781605352E2L,
2.2981886733594175366172E1L,
1.0000000000000000000000E0L,
];
// tan(PI/8)
enum real TAN_PI_8 = 4.1421356237309504880169e-1L;
// tan(3 * PI/8)
enum real TAN3_PI_8 = 2.41421356237309504880169L;
// Special cases.
if (x == 0.0)
return x;
if (isInfinity(x))
return copysign(PI_2, x);
// Make argument positive but save the sign.
bool sign = false;
if (signbit(x))
{
sign = true;
x = -x;
}
// Range reduction.
real y;
if (x > TAN3_PI_8)
{
y = PI_2;
x = -(1.0 / x);
}
else if (x > TAN_PI_8)
{
y = PI_4;
x = (x - 1.0)/(x + 1.0);
}
else
y = 0.0;
// Rational form in x^^2.
real z = x * x;
y = y + (poly(z, P) / poly(z, Q)) * z * x + x;
return (sign) ? -y : y;
}
}
/// ditto
double atan(double x) @safe pure nothrow @nogc { return atan(cast(real)x); }
/// ditto
float atan(float x) @safe pure nothrow @nogc { return atan(cast(real)x); }
unittest
{
assert(equalsDigit(atan(std.math.sqrt(3.0)), PI / 3, useDigits));
}
/***************
* Calculates the arc tangent of y / x,
* returning a value ranging from -$(PI) to $(PI).
*
* $(TABLE_SV
* $(TR $(TH y) $(TH x) $(TH atan(y, x)))
* $(TR $(TD $(NAN)) $(TD anything) $(TD $(NAN)) )
* $(TR $(TD anything) $(TD $(NAN)) $(TD $(NAN)) )
* $(TR $(TD $(PLUSMN)0.0) $(TD $(GT)0.0) $(TD $(PLUSMN)0.0) )
* $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) $(TD $(PLUSMN)0.0) )
* $(TR $(TD $(PLUSMN)0.0) $(TD $(LT)0.0) $(TD $(PLUSMN)$(PI)))
* $(TR $(TD $(PLUSMN)0.0) $(TD -0.0) $(TD $(PLUSMN)$(PI)))
* $(TR $(TD $(GT)0.0) $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
* $(TR $(TD $(LT)0.0) $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) )
* $(TR $(TD $(GT)0.0) $(TD $(INFIN)) $(TD $(PLUSMN)0.0) )
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything) $(TD $(PLUSMN)$(PI)/2))
* $(TR $(TD $(GT)0.0) $(TD -$(INFIN)) $(TD $(PLUSMN)$(PI)) )
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN)) $(TD $(PLUSMN)$(PI)/4))
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN)) $(TD $(PLUSMN)3$(PI)/4))
* )
*/
real atan2(real y, real x) @trusted pure nothrow @nogc
{
version(InlineAsm_X86_Any)
{
version (Win64)
{
asm pure nothrow @nogc {
naked;
fld real ptr [RDX]; // y
fld real ptr [RCX]; // x
fpatan;
ret;
}
}
else
{
asm pure nothrow @nogc {
fld y;
fld x;
fpatan;
}
}
}
else
{
// Special cases.
if (isNaN(x) || isNaN(y))
return real.nan;
if (y == 0.0)
{
if (x >= 0 && !signbit(x))
return copysign(0, y);
else
return copysign(PI, y);
}
if (x == 0.0)
return copysign(PI_2, y);
if (isInfinity(x))
{
if (signbit(x))
{
if (isInfinity(y))
return copysign(3*PI_4, y);
else
return copysign(PI, y);
}
else
{
if (isInfinity(y))
return copysign(PI_4, y);
else
return copysign(0.0, y);
}
}
if (isInfinity(y))
return copysign(PI_2, y);
// Call atan and determine the quadrant.
real z = atan(y / x);
if (signbit(x))
{
if (signbit(y))
z = z - PI;
else
z = z + PI;
}
if (z == 0.0)
return copysign(z, y);
return z;
}
}
/// ditto
double atan2(double y, double x) @safe pure nothrow @nogc
{
return atan2(cast(real)y, cast(real)x);
}
/// ditto
float atan2(float y, float x) @safe pure nothrow @nogc
{
return atan2(cast(real)y, cast(real)x);
}
unittest
{
assert(equalsDigit(atan2(1.0L, std.math.sqrt(3.0L)), PI / 6, useDigits));
}
/***********************************
* Calculates the hyperbolic cosine of x.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH cosh(x)) $(TH invalid?))
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
* )
*/
real cosh(real x) @safe pure nothrow @nogc
{
// cosh = (exp(x)+exp(-x))/2.
// The naive implementation works correctly.
real y = exp(x);
return (y + 1.0/y) * 0.5;
}
/// ditto
double cosh(double x) @safe pure nothrow @nogc { return cosh(cast(real)x); }
/// ditto
float cosh(float x) @safe pure nothrow @nogc { return cosh(cast(real)x); }
unittest
{
assert(equalsDigit(cosh(1.0), (E + 1.0 / E) / 2, useDigits));
}
/***********************************
* Calculates the hyperbolic sine of x.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH sinh(x)) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
* )
*/
real sinh(real x) @safe pure nothrow @nogc
{
// 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)
{
return copysign(0.5 * exp(fabs(x)), x);
}
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); }
unittest
{
assert(equalsDigit(sinh(1.0), (E - 1.0 / E) / 2, useDigits));
}
/***********************************
* Calculates the hyperbolic tangent of x.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH tanh(x)) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) )
* $(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);
}
real y = expm1(2*x);
return y / (y + 2);
}
/// ditto
double tanh(double x) @safe pure nothrow @nogc { return tanh(cast(real)x); }
/// ditto
float tanh(float x) @safe pure nothrow @nogc { return tanh(cast(real)x); }
unittest
{
assert(equalsDigit(tanh(1.0), sinh(1.0) / cosh(1.0), 15));
}
package:
/* Returns cosh(x) + I * sinh(x)
* Only one call to exp() is performed.
*/
creal coshisinh(real x) @safe pure nothrow @nogc
{
// See comments for cosh, sinh.
if (fabs(x) > real.mant_dig * LN2)
{
real y = exp(fabs(x));
return y * 0.5 + 0.5i * copysign(y, x);
}
else
{
real y = expm1(x);
return (y + 1.0 + 1.0/(y + 1.0)) * 0.5 + 0.5i * y / (y+1) * (y+2);
}
}
@safe pure nothrow @nogc unittest
{
creal c = coshisinh(3.0L);
assert(c.re == cosh(3.0L));
assert(c.im == sinh(3.0L));
}
public:
/***********************************
* Calculates the inverse hyperbolic cosine of x.
*
* Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
*
* $(TABLE_DOMRG
* $(DOMAIN 1..$(INFIN))
* $(RANGE 1..log(real.max), $(INFIN)) )
* $(TABLE_SV
* $(SVH x, acosh(x) )
* $(SV $(NAN), $(NAN) )
* $(SV $(LT)1, $(NAN) )
* $(SV 1, 0 )
* $(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));
}
/// ditto
double acosh(double x) @safe pure nothrow @nogc { return acosh(cast(real)x); }
/// ditto
float acosh(float x) @safe pure nothrow @nogc { return acosh(cast(real)x); }
unittest
{
assert(isNaN(acosh(0.9)));
assert(isNaN(acosh(real.nan)));
assert(acosh(1.0)==0.0);
assert(acosh(real.infinity) == real.infinity);
assert(isNaN(acosh(0.5)));
assert(equalsDigit(acosh(cosh(3.0)), 3, useDigits));
}
/***********************************
* Calculates the inverse hyperbolic sine of x.
*
* Mathematically,
* ---------------
* asinh(x) = log( x + sqrt( x*x + 1 )) // if x >= +0
* asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
* -------------
*
* $(TABLE_SV
* $(SVH x, asinh(x) )
* $(SV $(NAN), $(NAN) )
* $(SV $(PLUSMN)0, $(PLUSMN)0 )
* $(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);
}
/// ditto
double asinh(double x) @safe pure nothrow @nogc { return asinh(cast(real)x); }
/// ditto
float asinh(float x) @safe pure nothrow @nogc { return asinh(cast(real)x); }
unittest
{
assert(isIdentical(asinh(0.0), 0.0));
assert(isIdentical(asinh(-0.0), -0.0));
assert(asinh(real.infinity) == real.infinity);
assert(asinh(-real.infinity) == -real.infinity);
assert(isNaN(asinh(real.nan)));
assert(equalsDigit(asinh(sinh(3.0)), 3, useDigits));
}
/***********************************
* Calculates the inverse hyperbolic tangent of x,
* returning a value from ranging from -1 to 1.
*
* Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
*
*
* $(TABLE_DOMRG
* $(DOMAIN -$(INFIN)..$(INFIN))
* $(RANGE -1..1) )
* $(TABLE_SV
* $(SVH x, acosh(x) )
* $(SV $(NAN), $(NAN) )
* $(SV $(PLUSMN)0, $(PLUSMN)0)
* $(SV -$(INFIN), -0)
* )
*/
real atanh(real x) @safe pure nothrow @nogc
{
// log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
return 0.5 * log1p( 2 * x / (1 - x) );
}
/// ditto
double atanh(double x) @safe pure nothrow @nogc { return atanh(cast(real)x); }
/// ditto
float atanh(float x) @safe pure nothrow @nogc { return atanh(cast(real)x); }
unittest
{
assert(isIdentical(atanh(0.0), 0.0));
assert(isIdentical(atanh(-0.0),-0.0));
assert(isNaN(atanh(real.nan)));
assert(isNaN(atanh(-real.infinity)));
assert(atanh(0.0) == 0);
assert(equalsDigit(atanh(tanh(0.5L)), 0.5, useDigits));
}
/*****************************************
* Returns x rounded to a long value using the current rounding mode.
* If the integer value of x is
* greater than long.max, the result is
* indeterminate.
*/
long rndtol(real x) @nogc @safe pure nothrow; /* intrinsic */
//FIXME
///ditto
long rndtol(double x) @safe pure nothrow @nogc { return rndtol(cast(real)x); }
//FIXME
///ditto
long rndtol(float x) @safe pure nothrow @nogc { return rndtol(cast(real)x); }
/*****************************************
* Returns x rounded to a long value using the FE_TONEAREST rounding mode.
* If the integer value of x is
* greater than long.max, the result is
* indeterminate.
*/
extern (C) real rndtonl(real x);
/***************************************
* Compute square root of x.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH sqrt(x)) $(TH invalid?))
* $(TR $(TD -0.0) $(TD -0.0) $(TD no))
* $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD yes))
* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
* )
*/
float sqrt(float x) @nogc @safe pure nothrow; /* intrinsic */
/// ditto
double sqrt(double x) @nogc @safe pure nothrow; /* intrinsic */
/// ditto
real sqrt(real x) @nogc @safe pure nothrow; /* intrinsic */
@safe pure nothrow @nogc unittest
{
//ctfe
enum ZX80 = sqrt(7.0f);
enum ZX81 = sqrt(7.0);
enum ZX82 = sqrt(7.0L);
assert(isNaN(sqrt(-1.0f)));
assert(isNaN(sqrt(-1.0)));
assert(isNaN(sqrt(-1.0L)));
}
creal sqrt(creal z) @nogc @safe pure nothrow
{
creal c;
real x,y,w,r;
if (z == 0)
{
c = 0 + 0i;
}
else
{
real z_re = z.re;
real z_im = z.im;
x = fabs(z_re);
y = fabs(z_im);
if (x >= y)
{
r = y / x;
w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
}
else
{
r = x / y;
w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
}
if (z_re >= 0)
{
c = w + (z_im / (w + w)) * 1.0i;
}
else
{
if (z_im < 0)
w = -w;
c = z_im / (w + w) + w * 1.0i;
}
}
return c;
}
/**
* Calculates e$(SUPERSCRIPT x).
*
* $(TABLE_SV
* $(TR $(TH x) $(TH e$(SUPERSCRIPT x)) )
* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
* $(TR $(TD -$(INFIN)) $(TD +0.0) )
* $(TR $(TD $(NAN)) $(TD $(NAN)) )
* )
*/
real exp(real x) @trusted pure nothrow @nogc
{
version(D_InlineAsm_X86)
{
// e^^x = 2^^(LOG2E*x)
// (This is valid because the overflow & underflow limits for exp
// and exp2 are so similar).
return exp2(LOG2E*x);
}
else version(D_InlineAsm_X86_64)
{
// e^^x = 2^^(LOG2E*x)
// (This is valid because the overflow & underflow limits for exp
// and exp2 are so similar).
return exp2(LOG2E*x);
}
else
{
// Coefficients for exp(x)
static immutable real[3] P = [
9.9999999999999999991025E-1L,
3.0299440770744196129956E-2L,
1.2617719307481059087798E-4L,
];
static immutable real[4] Q = [
2.0000000000000000000897E0L,
2.2726554820815502876593E-1L,
2.5244834034968410419224E-3L,
3.0019850513866445504159E-6L,
];
// C1 + C2 = LN2.
enum real C1 = 6.9314575195312500000000E-1L;
enum real C2 = 1.428606820309417232121458176568075500134E-6L;
// Overflow and Underflow limits.
static if (real.mant_dig == 64) // 80-bit reals
{
enum real OF = 1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384) ≈ 0x1.62e42fefa39efp+13
enum real UF = -1.1398805384308300613366E4L; // ln(2^-16445) ≈ -0x1.6436716d5406ep+13
}
else static if (real.mant_dig == 53) // 64-bit reals
{
enum real OF = 0x1.62e42fefa39efp+9L; // ln((1-2^-53) * 2^1024) = 709.78271...
enum real UF = -0x1.74385446d71c3p+9L; // ln(2^-1074) = -744.44007...
}
else
static assert(0, "No over-/underflow limits for real type!");
// Special cases.
// FIXME: set IEEE flags accordingly
if (isNaN(x))
return x;
if (x > OF)
return real.infinity;
if (x < UF)
return 0.0;
// Express: e^^x = e^^g * 2^^n
// = e^^g * e^^(n * LOG2E)
// = e^^(g + n * LOG2E)
int n = cast(int)floor(LOG2E * x + 0.5);
x -= n * C1;
x -= n * C2;
// Rational approximation for exponential of the fractional part:
// e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
real xx = x * x;
real px = x * poly(xx, P);
x = px / (poly(xx, Q) - px);
x = 1.0 + ldexp(x, 1);
// Scale by power of 2.
x = ldexp(x, n);
return x;
}
}
/// ditto
double exp(double x) @safe pure nothrow @nogc { return exp(cast(real)x); }
/// ditto
float exp(float x) @safe pure nothrow @nogc { return exp(cast(real)x); }
unittest
{
assert(equalsDigit(exp(3.0L), E * E * E, useDigits));
}
/**
* Calculates the value of the natural logarithm base (e)
* raised to the power of x, minus 1.
*
* For very small x, expm1(x) is more accurate
* than exp(x)-1.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH e$(SUPERSCRIPT x)-1) )
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
* $(TR $(TD -$(INFIN)) $(TD -1.0) )
* $(TR $(TD $(NAN)) $(TD $(NAN)) )
* )
*/
real expm1(real x) @trusted pure nothrow @nogc
{
version(D_InlineAsm_X86)
{
enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
asm pure nothrow @nogc
{
/* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
* Author: Don Clugston.
*
* expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
* = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
* and 2ym1 = (2^^(y-rndint(y))-1).
* If 2rndy < 0.5*real.epsilon, result is -1.
* Implementation is otherwise the same as for exp2()
*/
naked;
fld real ptr [ESP+4] ; // x
mov AX, [ESP+4+8]; // AX = exponent and sign
sub ESP, 12+8; // Create scratch space on the stack
// [ESP,ESP+2] = scratchint
// [ESP+4..+6, +8..+10, +10] = scratchreal
// set scratchreal mantissa = 1.0
mov dword ptr [ESP+8], 0;
mov dword ptr [ESP+8+4], 0x80000000;
and AX, 0x7FFF; // drop sign bit
cmp AX, 0x401D; // avoid InvalidException in fist
jae L_extreme;
fldl2e;
fmulp ST(1), ST; // y = x*log2(e)
fist dword ptr [ESP]; // scratchint = rndint(y)
fisub dword ptr [ESP]; // y - rndint(y)
// and now set scratchreal exponent
mov EAX, [ESP];
add EAX, 0x3fff;
jle short L_largenegative;
cmp EAX,0x8000;
jge short L_largepositive;
mov [ESP+8+8],AX;
f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
fmul ST(1), ST; // ST=2rndy, ST(1)=2rndy*2ym1
fld1;
fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
add ESP,12+8;
ret PARAMSIZE;
L_extreme: // Extreme exponent. X is very large positive, very
// large negative, infinity, or NaN.
fxam;
fstsw AX;
test AX, 0x0400; // NaN_or_zero, but we already know x!=0
jz L_was_nan; // if x is NaN, returns x
test AX, 0x0200;
jnz L_largenegative;
L_largepositive:
// Set scratchreal = real.max.
// squaring it will create infinity, and set overflow flag.
mov word ptr [ESP+8+8], 0x7FFE;
fstp ST(0);
fld real ptr [ESP+8]; // load scratchreal
fmul ST(0), ST; // square it, to create havoc!
L_was_nan:
add ESP,12+8;
ret PARAMSIZE;
L_largenegative:
fstp ST(0);
fld1;
fchs; // return -1. Underflow flag is not set.
add ESP,12+8;
ret PARAMSIZE;
}
}
else version(D_InlineAsm_X86_64)
{
asm pure nothrow @nogc
{
naked;
}
version (Win64)
{
asm pure nothrow @nogc
{
fld real ptr [RCX]; // x
mov AX,[RCX+8]; // AX = exponent and sign
}
}
else
{
asm pure nothrow @nogc
{
fld real ptr [RSP+8]; // x
mov AX,[RSP+8+8]; // AX = exponent and sign
}
}
asm pure nothrow @nogc
{
/* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
* Author: Don Clugston.
*
* expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
* = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
* and 2ym1 = (2^(y-rndint(y))-1).
* If 2rndy < 0.5*real.epsilon, result is -1.
* Implementation is otherwise the same as for exp2()
*/
sub RSP, 24; // Create scratch space on the stack
// [RSP,RSP+2] = scratchint
// [RSP+4..+6, +8..+10, +10] = scratchreal
// set scratchreal mantissa = 1.0
mov dword ptr [RSP+8], 0;
mov dword ptr [RSP+8+4], 0x80000000;
and AX, 0x7FFF; // drop sign bit
cmp AX, 0x401D; // avoid InvalidException in fist
jae L_extreme;
fldl2e;
fmul ; // y = x*log2(e)
fist dword ptr [RSP]; // scratchint = rndint(y)
fisub dword ptr [RSP]; // y - rndint(y)
// and now set scratchreal exponent
mov EAX, [RSP];
add EAX, 0x3fff;
jle short L_largenegative;
cmp EAX,0x8000;
jge short L_largepositive;
mov [RSP+8+8],AX;
f2xm1; // 2^(y-rndint(y)) -1
fld real ptr [RSP+8] ; // 2^rndint(y)
fmul ST(1), ST;
fld1;
fsubp ST(1), ST;
fadd;
add RSP,24;
ret;
L_extreme: // Extreme exponent. X is very large positive, very
// large negative, infinity, or NaN.
fxam;
fstsw AX;
test AX, 0x0400; // NaN_or_zero, but we already know x!=0
jz L_was_nan; // if x is NaN, returns x
test AX, 0x0200;
jnz L_largenegative;
L_largepositive:
// Set scratchreal = real.max.
// squaring it will create infinity, and set overflow flag.
mov word ptr [RSP+8+8], 0x7FFE;
fstp ST(0);
fld real ptr [RSP+8]; // load scratchreal
fmul ST(0), ST; // square it, to create havoc!
L_was_nan:
add RSP,24;
ret;
L_largenegative:
fstp ST(0);
fld1;
fchs; // return -1. Underflow flag is not set.
add RSP,24;
ret;
}
}
else
{
// Coefficients for exp(x) - 1
static immutable real[5] P = [
-1.586135578666346600772998894928250240826E4L,
2.642771505685952966904660652518429479531E3L,
-3.423199068835684263987132888286791620673E2L,
1.800826371455042224581246202420972737840E1L,
-5.238523121205561042771939008061958820811E-1L,
];
static immutable real[6] Q = [
-9.516813471998079611319047060563358064497E4L,
3.964866271411091674556850458227710004570E4L,
-7.207678383830091850230366618190187434796E3L,
7.206038318724600171970199625081491823079E2L,
-4.002027679107076077238836622982900945173E1L,
1.000000000000000000000000000000000000000E0L,
];
// C1 + C2 = LN2.
enum real C1 = 6.9314575195312500000000E-1L;
enum real C2 = 1.4286068203094172321215E-6L;
// Overflow and Underflow limits.
enum real OF = 1.1356523406294143949492E4L;
enum real UF = -4.5054566736396445112120088E1L;
// Special cases.
if (x > OF)
return real.infinity;
if (x == 0.0)
return x;
if (x < UF)
return -1.0;
// Express x = LN2 (n + remainder), remainder not exceeding 1/2.
int n = cast(int)floor(0.5 + x / LN2);
x -= n * C1;
x -= n * C2;
// Rational approximation:
// exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x)
real px = x * poly(x, P);
real qx = poly(x, Q);
real xx = x * x;
qx = x + (0.5 * xx + xx * px / qx);
// We have qx = exp(remainder LN2) - 1, so:
// exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1.
px = ldexp(1.0, n);
x = px * qx + (px - 1.0);
return x;
}
}
/**
* Calculates 2$(SUPERSCRIPT x).
*
* $(TABLE_SV
* $(TR $(TH x) $(TH exp2(x)) )
* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
* $(TR $(TD -$(INFIN)) $(TD +0.0) )
* $(TR $(TD $(NAN)) $(TD $(NAN)) )
* )
*/
real exp2(real x) @nogc @trusted pure nothrow
{
version(D_InlineAsm_X86)
{
enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
asm pure nothrow @nogc
{
/* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
* Author: Don Clugston.
*
* exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
* The trick for high performance is to avoid the fscale(28cycles on core2),
* frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
*
* We can do frndint by using fist. BUT we can't use it for huge numbers,
* because it will set the Invalid Operation flag if overflow or NaN occurs.
* Fortunately, whenever this happens the result would be zero or infinity.
*
* We can perform fscale by directly poking into the exponent. BUT this doesn't
* work for the (very rare) cases where the result is subnormal. So we fall back
* to the slow method in that case.
*/
naked;
fld real ptr [ESP+4] ; // x
mov AX, [ESP+4+8]; // AX = exponent and sign
sub ESP, 12+8; // Create scratch space on the stack
// [ESP,ESP+2] = scratchint
// [ESP+4..+6, +8..+10, +10] = scratchreal
// set scratchreal mantissa = 1.0
mov dword ptr [ESP+8], 0;
mov dword ptr [ESP+8+4], 0x80000000;
and AX, 0x7FFF; // drop sign bit
cmp AX, 0x401D; // avoid InvalidException in fist
jae L_extreme;
fist dword ptr [ESP]; // scratchint = rndint(x)
fisub dword ptr [ESP]; // x - rndint(x)
// and now set scratchreal exponent
mov EAX, [ESP];
add EAX, 0x3fff;
jle short L_subnormal;
cmp EAX,0x8000;
jge short L_overflow;
mov [ESP+8+8],AX;
L_normal:
f2xm1;
fld1;
faddp ST(1), ST; // 2^^(x-rndint(x))
fld real ptr [ESP+8] ; // 2^^rndint(x)
add ESP,12+8;
fmulp ST(1), ST;
ret PARAMSIZE;
L_subnormal:
// Result will be subnormal.
// In this rare case, the simple poking method doesn't work.
// The speed doesn't matter, so use the slow fscale method.
fild dword ptr [ESP]; // scratchint
fld1;
fscale;
fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
fstp ST(0); // drop scratchint
jmp L_normal;
L_extreme: // Extreme exponent. X is very large positive, very
// large negative, infinity, or NaN.
fxam;
fstsw AX;
test AX, 0x0400; // NaN_or_zero, but we already know x!=0
jz L_was_nan; // if x is NaN, returns x
// set scratchreal = real.min_normal
// squaring it will return 0, setting underflow flag
mov word ptr [ESP+8+8], 1;
test AX, 0x0200;
jnz L_waslargenegative;
L_overflow:
// Set scratchreal = real.max.
// squaring it will create infinity, and set overflow flag.
mov word ptr [ESP+8+8], 0x7FFE;
L_waslargenegative:
fstp ST(0);
fld real ptr [ESP+8]; // load scratchreal
fmul ST(0), ST; // square it, to create havoc!
L_was_nan:
add ESP,12+8;
ret PARAMSIZE;
}
}
else version(D_InlineAsm_X86_64)
{
asm pure nothrow @nogc
{
naked;
}
version (Win64)
{
asm pure nothrow @nogc
{
fld real ptr [RCX]; // x
mov AX,[RCX+8]; // AX = exponent and sign
}
}
else
{
asm pure nothrow @nogc
{
fld real ptr [RSP+8]; // x
mov AX,[RSP+8+8]; // AX = exponent and sign
}
}
asm pure nothrow @nogc
{
/* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
* Author: Don Clugston.
*
* exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
* The trick for high performance is to avoid the fscale(28cycles on core2),
* frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
*
* We can do frndint by using fist. BUT we can't use it for huge numbers,
* because it will set the Invalid Operation flag is overflow or NaN occurs.
* Fortunately, whenever this happens the result would be zero or infinity.
*
* We can perform fscale by directly poking into the exponent. BUT this doesn't
* work for the (very rare) cases where the result is subnormal. So we fall back
* to the slow method in that case.
*/
sub RSP, 24; // Create scratch space on the stack
// [RSP,RSP+2] = scratchint
// [RSP+4..+6, +8..+10, +10] = scratchreal
// set scratchreal mantissa = 1.0
mov dword ptr [RSP+8], 0;
mov dword ptr [RSP+8+4], 0x80000000;
and AX, 0x7FFF; // drop sign bit
cmp AX, 0x401D; // avoid InvalidException in fist
jae L_extreme;
fist dword ptr [RSP]; // scratchint = rndint(x)
fisub dword ptr [RSP]; // x - rndint(x)
// and now set scratchreal exponent
mov EAX, [RSP];
add EAX, 0x3fff;
jle short L_subnormal;
cmp EAX,0x8000;
jge short L_overflow;
mov [RSP+8+8],AX;
L_normal:
f2xm1;
fld1;
fadd; // 2^(x-rndint(x))
fld real ptr [RSP+8] ; // 2^rndint(x)
add RSP,24;
fmulp ST(1), ST;
ret;
L_subnormal:
// Result will be subnormal.
// In this rare case, the simple poking method doesn't work.
// The speed doesn't matter, so use the slow fscale method.
fild dword ptr [RSP]; // scratchint
fld1;
fscale;
fstp real ptr [RSP+8]; // scratchreal = 2^scratchint
fstp ST(0); // drop scratchint
jmp L_normal;
L_extreme: // Extreme exponent. X is very large positive, very
// large negative, infinity, or NaN.
fxam;
fstsw AX;
test AX, 0x0400; // NaN_or_zero, but we already know x!=0
jz L_was_nan; // if x is NaN, returns x
// set scratchreal = real.min
// squaring it will return 0, setting underflow flag
mov word ptr [RSP+8+8], 1;
test AX, 0x0200;
jnz L_waslargenegative;
L_overflow:
// Set scratchreal = real.max.
// squaring it will create infinity, and set overflow flag.
mov word ptr [RSP+8+8], 0x7FFE;
L_waslargenegative:
fstp ST(0);
fld real ptr [RSP+8]; // load scratchreal
fmul ST(0), ST; // square it, to create havoc!
L_was_nan:
add RSP,24;
ret;
}
}
else
{
// Coefficients for exp2(x)
static immutable real[3] P = [
2.0803843631901852422887E6L,
3.0286971917562792508623E4L,
6.0614853552242266094567E1L,
];
static immutable real[4] Q = [
6.0027204078348487957118E6L,
3.2772515434906797273099E5L,
1.7492876999891839021063E3L,
1.0000000000000000000000E0L,
];
// Overflow and Underflow limits.
enum real OF = 16384.0L;
enum real UF = -16382.0L;
// Special cases.
if (isNaN(x))
return x;
if (x > OF)
return real.infinity;
if (x < UF)
return 0.0;
// Separate into integer and fractional parts.
int n = cast(int)floor(x + 0.5);
x -= n;
// Rational approximation:
// exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
real xx = x * x;
real px = x * poly(xx, P);
x = px / (poly(xx, Q) - px);
x = 1.0 + ldexp(x, 1);
// Scale by power of 2.
x = ldexp(x, n);
return x;
}
}
///
unittest
{
assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1);
assert(exp2(8.0L) == 256.0);
assert(exp2(-9.0L)== 1.0L/512.0);
}
unittest
{
version(CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented
{
assert( core.stdc.math.exp2f(0.0f) == 1 );
assert( core.stdc.math.exp2 (0.0) == 1 );
assert( core.stdc.math.exp2l(0.0L) == 1 );
}
}
unittest
{
FloatingPointControl ctrl;
if(FloatingPointControl.hasExceptionTraps)
ctrl.disableExceptions(FloatingPointControl.allExceptions);
ctrl.rounding = FloatingPointControl.roundToNearest;
static if (real.mant_dig == 64) // 80-bit reals
{
static immutable real[2][] exptestpoints =
[ // x exp(x)
[ 1.0L, E ],
[ 0.5L, 0x1.a61298e1e069bc97p+0L ],
[ 3.0L, E*E*E ],
[ 0x1.1p+13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow
[ 0x1.7p+13L, real.infinity ], // close overflow
[ 0x1p+80L, real.infinity ], // far overflow
[ real.infinity, real.infinity ],
[-0x1.18p+13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow
[-0x1.625p+13L, 0x1.a6bd68a39d11f35cp-16358L ], // ditto
[-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L ], // near underflow - subnormal
[-0x1.643p+13L, 0x1p-16444L ], // ditto
[-0x1.645p+13L, 0 ], // close underflow
[-0x1p+30L, 0 ], // far underflow
];
}
else static if (real.mant_dig == 53) // 64-bit reals
{
static immutable real[2][] exptestpoints =
[ // x, exp(x)
[ 1.0L, E ],
[ 0.5L, 0x1.a61298e1e069cp+0L ],
[ 3.0L, E*E*E ],
[ 0x1.6p+9L, 0x1.93bf4ec282efbp+1015L ], // near overflow
[ 0x1.7p+9L, real.infinity ], // close overflow
[ 0x1p+80L, real.infinity ], // far overflow
[ real.infinity, real.infinity ],
[-0x1.6p+9L, 0x1.44a3824e5285fp-1016L ], // near underflow
[-0x1.64p+9L, 0x0.06f84920bb2d3p-1022L ], // near underflow - subnormal
[-0x1.743p+9L, 0x0.0000000000001p-1022L ], // ditto
[-0x1.8p+9L, 0 ], // close underflow
[-0x1p30L, 0 ], // far underflow
];
}
else
static assert(0, "No exp() tests for real type!");
const minEqualMantissaBits = real.mant_dig - 2;
real x;
IeeeFlags f;
foreach (ref pair; exptestpoints)
{
resetIeeeFlags();
x = exp(pair[0]);
f = ieeeFlags;
assert(feqrel(x, pair[1]) >= minEqualMantissaBits);
version (IeeeFlagsSupport)
{
// Check the overflow bit
if (x == real.infinity)
{
// don't care about the overflow bit if input was inf
// (e.g., the LLVM intrinsic doesn't set it on Linux x86_64)
assert(pair[0] == real.infinity || f.overflow);
}
else
assert(!f.overflow);
// Check the underflow bit
assert(f.underflow == (fabs(x) < real.min_normal));
// Invalid and div by zero shouldn't be affected.
assert(!f.invalid);
assert(!f.divByZero);
}
}
// Ideally, exp(0) would not set the inexact flag.
// Unfortunately, fldl2e sets it!
// So it's not realistic to avoid setting it.
assert(exp(0.0L) == 1.0);
// NaN propagation. Doesn't set flags, bcos was already NaN.
resetIeeeFlags();
x = exp(real.nan);
f = ieeeFlags;
assert(isIdentical(abs(x), real.nan));
assert(f.flags == 0);
resetIeeeFlags();
x = exp(-real.nan);
f = ieeeFlags;
assert(isIdentical(abs(x), real.nan));
assert(f.flags == 0);
x = exp(NaN(0x123));
assert(isIdentical(x, NaN(0x123)));
// High resolution test
assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6D_33Fp+0L);
}
/**
* Calculate cos(y) + i sin(y).
*
* On many CPUs (such as x86), this is a very efficient operation;
* almost twice as fast as calculating sin(y) and cos(y) separately,
* and is the preferred method when both are required.
*/
creal expi(real y) @trusted pure nothrow @nogc
{
version(InlineAsm_X86_Any)
{
version (Win64)
{
asm pure nothrow @nogc
{
naked;
fld real ptr [ECX];
fsincos;
fxch ST(1), ST(0);
ret;
}
}
else
{
asm pure nothrow @nogc
{
fld y;
fsincos;
fxch ST(1), ST(0);
}
}
}
else
{
return cos(y) + sin(y)*1i;
}
}
///
@safe pure nothrow @nogc unittest
{
assert(expi(1.3e5L) == cos(1.3e5L) + sin(1.3e5L) * 1i);
assert(expi(0.0L) == 1L + 0.0Li);
}
/*********************************************************************
* Separate floating point value into significand and exponent.
*
* Returns:
* Calculate and return $(I x) and $(I exp) such that
* value =$(I x)*2$(SUPERSCRIPT exp) and
* .5 $(LT)= |$(I x)| $(LT) 1.0
*
* $(I x) has same sign as value.
*
* $(TABLE_SV
* $(TR $(TH value) $(TH returns) $(TH exp))
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0))
* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max))
* $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min))
* $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
* )
*/
T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc
if(isFloatingPoint!T)
{
Unqual!T vf = value;
ushort* vu = cast(ushort*)&vf;
static if(is(Unqual!T == float))
int* vi = cast(int*)&vf;
else
long* vl = cast(long*)&vf;
int ex;
alias F = floatTraits!T;
ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
static if (F.realFormat == RealFormat.ieeeExtended)
{
if (ex)
{ // If exponent is non-zero
if (ex == F.EXPMASK) // infinity or NaN
{
if (*vl & 0x7FFF_FFFF_FFFF_FFFF) // NaN
{
*vl |= 0xC000_0000_0000_0000; // convert NaNS to NaNQ
exp = int.min;
}
else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity
exp = int.min;
else // positive infinity
exp = int.max;
}
else
{
exp = ex - F.EXPBIAS;
vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE;
}
}
else if (!*vl)
{
// vf is +-0.0
exp = 0;
}
else
{
// subnormal
vf *= F.RECIP_EPSILON;
ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
exp = ex - F.EXPBIAS - T.mant_dig + 1;
vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE;
}
return vf;
}
else static if (F.realFormat == RealFormat.ieeeQuadruple)
{
if (ex) // If exponent is non-zero
{
if (ex == F.EXPMASK)
{ // infinity or NaN
if (vl[MANTISSA_LSB] |
( vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN
{
// convert NaNS to NaNQ
vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000;
exp = int.min;
}
else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity
exp = int.min;
else // positive infinity
exp = int.max;
}
else
{
exp = ex - F.EXPBIAS;
vu[F.EXPPOS_SHORT] =
cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE);
}
}
else if ((vl[MANTISSA_LSB]
|(vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
{
// vf is +-0.0
exp = 0;
}
else
{
// subnormal
vf *= F.RECIP_EPSILON;
ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
exp = ex - F.EXPBIAS - T.mant_dig + 1;
vu[F.EXPPOS_SHORT] =
cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE);
}
return vf;
}
else static if (F.realFormat == RealFormat.ieeeDouble)
{
if (ex) // If exponent is non-zero
{
if (ex == F.EXPMASK) // infinity or NaN
{
if (*vl == 0x7FF0_0000_0000_0000) // positive infinity
{
exp = int.max;
}
else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity
exp = int.min;
else
{ // NaN
*vl |= 0x0008_0000_0000_0000; // convert NaNS to NaNQ
exp = int.min;
}
}
else
{
exp = (ex - F.EXPBIAS) >> 4;
vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0);
}
}
else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF))
{
// vf is +-0.0
exp = 0;
}
else
{
// subnormal
vf *= F.RECIP_EPSILON;
ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1;
vu[F.EXPPOS_SHORT] =
cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FE0);
}
return vf;
}
else static if (F.realFormat == RealFormat.ieeeSingle)
{
if (ex) // If exponent is non-zero
{
if (ex == F.EXPMASK) // infinity or NaN
{
if (*vi == 0x7F80_0000) // positive infinity
{
exp = int.max;
}
else if (*vi == 0xFF80_0000) // negative infinity
exp = int.min;
else
{ // NaN
*vi |= 0x0040_0000; // convert NaNS to NaNQ
exp = int.min;
}
}
else
{
exp = (ex - F.EXPBIAS) >> 7;
vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00);
}
}
else if (!(*vi & 0x7FFF_FFFF))
{
// vf is +-0.0
exp = 0;
}
else
{
// subnormal
vf *= F.RECIP_EPSILON;
ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1;
vu[F.EXPPOS_SHORT] =
cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3F00);
}
return vf;
}
else // static if (F.realFormat == RealFormat.ibmExtended)
{
assert (0, "frexp not implemented");
}
}
///
unittest
{
int exp;
real mantissa = frexp(123.456L, exp);
// check if values are equal to 19 decimal digits of precision
assert(equalsDigit(mantissa * pow(2.0L, cast(real)exp), 123.456L, 19));
assert(frexp(-real.nan, exp) && exp == int.min);
assert(frexp(real.nan, exp) && exp == int.min);
assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min);
assert(frexp(real.infinity, exp) == real.infinity && exp == int.max);
assert(frexp(-0.0, exp) == -0.0 && exp == 0);
assert(frexp(0.0, exp) == 0.0 && exp == 0);
}
unittest
{
import std.typetuple, std.typecons;
foreach (T; TypeTuple!(real, double, float))
{
Tuple!(T, T, int)[] vals = // x,frexp,exp
[
tuple(T(0.0), T( 0.0 ), 0),
tuple(T(-0.0), T( -0.0), 0),
tuple(T(1.0), T( .5 ), 1),
tuple(T(-1.0), T( -.5 ), 1),
tuple(T(2.0), T( .5 ), 2),
tuple(T(float.min_normal/2.0f), T(.5), -126),
tuple(T.infinity, T.infinity, int.max),
tuple(-T.infinity, -T.infinity, int.min),
tuple(T.nan, T.nan, int.min),
tuple(-T.nan, -T.nan, int.min),
];
foreach(elem; vals)
{
T x = elem[0];
T e = elem[1];
int exp = elem[2];
int eptr;
T v = frexp(x, eptr);
assert(isIdentical(e, v));
assert(exp == eptr);
}
static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
{
static T[3][] extendedvals = [ // x,frexp,exp
[0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal
[0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
[T.min_normal, .5, -16381],
[T.min_normal/2.0L, .5, -16382] // subnormal
];
foreach(elem; extendedvals)
{
T x = elem[0];
T e = elem[1];
int exp = cast(int)elem[2];
int eptr;
T v = frexp(x, eptr);
assert(isIdentical(e, v));
assert(exp == eptr);
}
}
}
}
unittest
{
import std.typetuple: TypeTuple;
void foo() {
foreach (T; TypeTuple!(real, double, float))
{
int exp;
const T a = 1;
immutable T b = 2;
auto c = frexp(a, exp);
auto d = frexp(b, exp);
}
}
}
/******************************************
* Extracts the exponent of x as a signed integral value.
*
* If x is not a special value, the result is the same as
* $(D cast(int)logb(x)).
*
* $(TABLE_SV
* $(TR $(TH x) $(TH ilogb(x)) $(TH Range error?))
* $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes))
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max) $(TD no))
* $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD no))
* )
*/
int ilogb(real x) @trusted nothrow @nogc
{
version (Win64_DMD_InlineAsm)
{
asm pure nothrow @nogc
{
naked ;
fld real ptr [RCX] ;
fxam ;
fstsw AX ;
and AH,0x45 ;
cmp AH,0x40 ;
jz Lzeronan ;
cmp AH,5 ;
jz Linfinity ;
cmp AH,1 ;
jz Lzeronan ;
fxtract ;
fstp ST(0) ;
fistp dword ptr 8[RSP] ;
mov EAX,8[RSP] ;
ret ;
Lzeronan:
mov EAX,0x80000000 ;
fstp ST(0) ;
ret ;
Linfinity:
mov EAX,0x7FFFFFFF ;
fstp ST(0) ;
ret ;
}
}
else version (CRuntime_Microsoft)
{
int res;
asm pure nothrow @nogc
{
fld real ptr [x] ;
fxam ;
fstsw AX ;
and AH,0x45 ;
cmp AH,0x40 ;
jz Lzeronan ;
cmp AH,5 ;
jz Linfinity ;
cmp AH,1 ;
jz Lzeronan ;
fxtract ;
fstp ST(0) ;
fistp res ;
mov EAX,res ;
jmp Ldone ;
Lzeronan:
mov EAX,0x80000000 ;
fstp ST(0) ;
jmp Ldone ;
Linfinity:
mov EAX,0x7FFFFFFF ;
fstp ST(0) ;
Ldone: ;
}
}
else
return core.stdc.math.ilogbl(x);
}
alias FP_ILOGB0 = core.stdc.math.FP_ILOGB0;
alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN;
@trusted nothrow @nogc unittest
{
assert(ilogb(real.nan) == FP_ILOGBNAN);
assert(ilogb(-real.nan) == FP_ILOGBNAN);
assert(ilogb(0.0) == FP_ILOGB0);
assert(ilogb(-0.0) == FP_ILOGB0);
assert(ilogb(real.infinity) == int.max);
assert(ilogb(-real.infinity) == int.max);
assert(ilogb(2.0) == 1);
assert(ilogb(2.0001) == 1);
assert(ilogb(1.9999) == 0);
assert(ilogb(0.5) == -1);
assert(ilogb(123.123) == 6);
assert(ilogb(-123.123) == 6);
assert(ilogb(0.123) == -4);
assert(ilogb(-double.min_normal) == -1022);
assert(ilogb(-float.min_normal) == -126);
// subnormals
assert(ilogb(nextUp(-double.min_normal)) == -1023);
assert(ilogb(nextUp(-0.0)) == -1074);
assert(ilogb(nextUp(-float.min_normal)) == -127);
assert(ilogb(nextUp(-0.0F)) == -149);
static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended) {
assert(ilogb(-real.min_normal) == -16382);
assert(ilogb(nextUp(-real.min_normal)) == -16383);
assert(ilogb(nextUp(-0.0L)) == -16445);
} else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) {
assert(ilogb(-real.min_normal) == -1022);
assert(ilogb(nextUp(-real.min_normal)) == -1023);
assert(ilogb(nextUp(-0.0L)) == -1074);
}
}
/*******************************************
* Compute n * 2$(SUPERSCRIPT exp)
* References: frexp
*/
real ldexp(real n, int exp) @nogc @safe pure nothrow; /* intrinsic */
//FIXME
///ditto
double ldexp(double n, int exp) @safe pure nothrow @nogc { return ldexp(cast(real)n, exp); }
//FIXME
///ditto
float ldexp(float n, int exp) @safe pure nothrow @nogc { return ldexp(cast(real)n, exp); }
///
@nogc @safe pure nothrow unittest
{
import std.typetuple;
foreach(T; TypeTuple!(float, double, real))
{
T r;
r = ldexp(3.0L, 3);
assert(r == 24);
r = ldexp(cast(T)3.0, cast(int) 3);
assert(r == 24);
T n = 3.0;
int exp = 3;
r = ldexp(n, exp);
assert(r == 24);
}
}
@safe pure nothrow @nogc unittest
{
static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
{
assert(ldexp(1.0L, -16384) == 0x1p-16384L);
assert(ldexp(1.0L, -16382) == 0x1p-16382L);
int x;
real n = frexp(0x1p-16384L, x);
assert(n==0.5L);
assert(x==-16383);
assert(ldexp(n, x)==0x1p-16384L);
}
else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
{
assert(ldexp(1.0L, -1024) == 0x1p-1024L);
assert(ldexp(1.0L, -1022) == 0x1p-1022L);
int x;
real n = frexp(0x1p-1024L, x);
assert(n==0.5L);
assert(x==-1023);
assert(ldexp(n, x)==0x1p-1024L);
}
else static assert(false, "Floating point type real not supported");
}
@safe pure nothrow @nogc unittest
{
assert(ldexp(1.0, -1024) == 0x1p-1024);
assert(ldexp(1.0, -1022) == 0x1p-1022);
int x;
double n = frexp(0x1p-1024, x);
assert(n==0.5);
assert(x==-1023);
assert(ldexp(n, x)==0x1p-1024);
}
@safe pure nothrow @nogc unittest
{
assert(ldexp(1.0f, -128) == 0x1p-128f);
assert(ldexp(1.0f, -126) == 0x1p-126f);
int x;
float n = frexp(0x1p-128f, x);
assert(n==0.5f);
assert(x==-127);
assert(ldexp(n, x)==0x1p-128f);
}
unittest
{
static real[3][] vals = // value,exp,ldexp
[
[ 0, 0, 0],
[ 1, 0, 1],
[ -1, 0, -1],
[ 1, 1, 2],
[ 123, 10, 125952],
[ real.max, int.max, real.infinity],
[ real.max, -int.max, 0],
[ real.min_normal, -int.max, 0],
];
int i;
for (i = 0; i < vals.length; i++)
{
real x = vals[i][0];
int exp = cast(int)vals[i][1];
real z = vals[i][2];
real l = ldexp(x, exp);
assert(equalsDigit(z, l, 7));
}
}
/**************************************
* Calculate the natural logarithm of x.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH log(x)) $(TH divide by 0?) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
* $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
* )
*/
real log(real x) @safe pure nothrow @nogc
{
version (INLINE_YL2X)
return yl2x(x, LN2);
else
{
// Coefficients for log(1 + x)
static immutable real[7] P = [
2.0039553499201281259648E1L,
5.7112963590585538103336E1L,
6.0949667980987787057556E1L,
2.9911919328553073277375E1L,
6.5787325942061044846969E0L,
4.9854102823193375972212E-1L,
4.5270000862445199635215E-5L,
];
static immutable real[7] Q = [
6.0118660497603843919306E1L,
2.1642788614495947685003E2L,
3.0909872225312059774938E2L,
2.2176239823732856465394E2L,
8.3047565967967209469434E1L,
1.5062909083469192043167E1L,
1.0000000000000000000000E0L,
];
// Coefficients for log(x)
static immutable real[4] R = [
-3.5717684488096787370998E1L,
1.0777257190312272158094E1L,
-7.1990767473014147232598E-1L,
1.9757429581415468984296E-3L,
];
static immutable real[4] S = [
-4.2861221385716144629696E2L,
1.9361891836232102174846E2L,
-2.6201045551331104417768E1L,
1.0000000000000000000000E0L,
];
// C1 + C2 = LN2.
enum real C1 = 6.9314575195312500000000E-1L;
enum real C2 = 1.4286068203094172321215E-6L;
// Special cases.
if (isNaN(x))
return x;
if (isInfinity(x) && !signbit(x))
return x;
if (x == 0.0)
return -real.infinity;
if (x < 0.0)
return real.nan;
// Separate mantissa from exponent.
// Note, frexp is used so that denormal numbers will be handled properly.
real y, z;
int exp;
x = frexp(x, exp);
// Logarithm using log(x) = z + z^^3 P(z) / Q(z),
// where z = 2(x - 1)/(x + 1)
if((exp > 2) || (exp < -2))
{
if(x < SQRT1_2)
{ // 2(2x - 1)/(2x + 1)
exp -= 1;
z = x - 0.5;
y = 0.5 * z + 0.5;
}
else
{ // 2(x - 1)/(x + 1)
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
z = x * x;
z = x * (z * poly(z, R) / poly(z, S));
z += exp * C2;
z += x;
z += exp * C1;
return z;
}
// Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
if (x < SQRT1_2)
{ // 2x - 1
exp -= 1;
x = ldexp(x, 1) - 1.0;
}
else
{
x = x - 1.0;
}
z = x * x;
y = x * (z * poly(x, P) / poly(x, Q));
y += exp * C2;
z = y - ldexp(z, -1);
// Note, the sum of above terms does not exceed x/4,
// so it contributes at most about 1/4 lsb to the error.
z += x;
z += exp * C1;
return z;
}
}
///
@safe pure nothrow @nogc unittest
{
assert(log(E) == 1);
}
/**************************************
* Calculate the base-10 logarithm of x.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH log10(x)) $(TH divide by 0?) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
* $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
* )
*/
real log10(real x) @safe pure nothrow @nogc
{
version (INLINE_YL2X)
return yl2x(x, LOG2);
else
{
// Coefficients for log(1 + x)
static immutable real[7] P = [
2.0039553499201281259648E1L,
5.7112963590585538103336E1L,
6.0949667980987787057556E1L,
2.9911919328553073277375E1L,
6.5787325942061044846969E0L,
4.9854102823193375972212E-1L,
4.5270000862445199635215E-5L,
];
static immutable real[7] Q = [
6.0118660497603843919306E1L,
2.1642788614495947685003E2L,
3.0909872225312059774938E2L,
2.2176239823732856465394E2L,
8.3047565967967209469434E1L,
1.5062909083469192043167E1L,
1.0000000000000000000000E0L,
];
// Coefficients for log(x)
static immutable real[4] R = [
-3.5717684488096787370998E1L,
1.0777257190312272158094E1L,
-7.1990767473014147232598E-1L,
1.9757429581415468984296E-3L,
];
static immutable real[4] S = [
-4.2861221385716144629696E2L,
1.9361891836232102174846E2L,
-2.6201045551331104417768E1L,
1.0000000000000000000000E0L,
];
// log10(2) split into two parts.
enum real L102A = 0.3125L;
enum real L102B = -1.14700043360188047862611052755069732318101185E-2L;
// log10(e) split into two parts.
enum real L10EA = 0.5L;
enum real L10EB = -6.570551809674817234887108108339491770560299E-2L;
// Special cases are the same as for log.
if (isNaN(x))
return x;
if (isInfinity(x) && !signbit(x))
return x;
if (x == 0.0)
return -real.infinity;
if (x < 0.0)
return real.nan;
// Separate mantissa from exponent.
// Note, frexp is used so that denormal numbers will be handled properly.
real y, z;
int exp;
x = frexp(x, exp);
// Logarithm using log(x) = z + z^^3 P(z) / Q(z),
// where z = 2(x - 1)/(x + 1)
if((exp > 2) || (exp < -2))
{
if(x < SQRT1_2)
{ // 2(2x - 1)/(2x + 1)
exp -= 1;
z = x - 0.5;
y = 0.5 * z + 0.5;
}
else
{ // 2(x - 1)/(x + 1)
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
z = x * x;
y = x * (z * poly(z, R) / poly(z, S));
goto Ldone;
}
// Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
if (x < SQRT1_2)
{ // 2x - 1
exp -= 1;
x = ldexp(x, 1) - 1.0;
}
else
x = x - 1.0;
z = x * x;
y = x * (z * poly(x, P) / poly(x, Q));
y = y - ldexp(z, -1);
// Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
// This sequence of operations is critical and it may be horribly
// defeated by some compiler optimizers.
Ldone:
z = y * L10EB;
z += x * L10EB;
z += exp * L102B;
z += y * L10EA;
z += x * L10EA;
z += exp * L102A;
return z;
}
}
///
@safe pure nothrow @nogc unittest
{
assert(fabs(log10(1000) - 3) < .000001);
}
/******************************************
* Calculates the natural logarithm of 1 + x.
*
* For very small x, log1p(x) will be more accurate than
* log(1 + x).
*
* $(TABLE_SV
* $(TR $(TH x) $(TH log1p(x)) $(TH divide by 0?) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) $(TD no))
* $(TR $(TD -1.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
* $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD no) $(TD yes))
* $(TR $(TD +$(INFIN)) $(TD -$(INFIN)) $(TD no) $(TD no))
* )
*/
real log1p(real x) @safe pure nothrow @nogc
{
version(INLINE_YL2X)
{
// On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5,
// ie if -0.29<=x<=0.414
return (fabs(x) <= 0.25) ? yl2xp1(x, LN2) : yl2x(x+1, LN2);
}
else
{
// Special cases.
if (isNaN(x) || x == 0.0)
return x;
if (isInfinity(x) && !signbit(x))
return x;
if (x == -1.0)
return -real.infinity;
if (x < -1.0)
return real.nan;
return log(x + 1.0);
}
}
/***************************************
* Calculates the base-2 logarithm of x:
* $(SUB log, 2)x
*
* $(TABLE_SV
* $(TR $(TH x) $(TH log2(x)) $(TH divide by 0?) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no) )
* $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes) )
* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no) )
* )
*/
real log2(real x) @safe pure nothrow @nogc
{
version (INLINE_YL2X)
return yl2x(x, 1);
else
{
// Coefficients for log(1 + x)
static immutable real[7] P = [
2.0039553499201281259648E1L,
5.7112963590585538103336E1L,
6.0949667980987787057556E1L,
2.9911919328553073277375E1L,
6.5787325942061044846969E0L,
4.9854102823193375972212E-1L,
4.5270000862445199635215E-5L,
];
static immutable real[7] Q = [
6.0118660497603843919306E1L,
2.1642788614495947685003E2L,
3.0909872225312059774938E2L,
2.2176239823732856465394E2L,
8.3047565967967209469434E1L,
1.5062909083469192043167E1L,
1.0000000000000000000000E0L,
];
// Coefficients for log(x)
static immutable real[4] R = [
-3.5717684488096787370998E1L,
1.0777257190312272158094E1L,
-7.1990767473014147232598E-1L,
1.9757429581415468984296E-3L,
];
static immutable real[4] S = [
-4.2861221385716144629696E2L,
1.9361891836232102174846E2L,
-2.6201045551331104417768E1L,
1.0000000000000000000000E0L,
];
// Special cases are the same as for log.
if (isNaN(x))
return x;
if (isInfinity(x) && !signbit(x))
return x;
if (x == 0.0)
return -real.infinity;
if (x < 0.0)
return real.nan;
// Separate mantissa from exponent.
// Note, frexp is used so that denormal numbers will be handled properly.
real y, z;
int exp;
x = frexp(x, exp);
// Logarithm using log(x) = z + z^^3 P(z) / Q(z),
// where z = 2(x - 1)/(x + 1)
if((exp > 2) || (exp < -2))
{
if(x < SQRT1_2)
{ // 2(2x - 1)/(2x + 1)
exp -= 1;
z = x - 0.5;
y = 0.5 * z + 0.5;
}
else
{ // 2(x - 1)/(x + 1)
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
z = x * x;
y = x * (z * poly(z, R) / poly(z, S));
goto Ldone;
}
// Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
if (x < SQRT1_2)
{ // 2x - 1
exp -= 1;
x = ldexp(x, 1) - 1.0;
}
else
x = x - 1.0;
z = x * x;
y = x * (z * poly(x, P) / poly(x, Q));
y = y - ldexp(z, -1);
// Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
// This sequence of operations is critical and it may be horribly
// defeated by some compiler optimizers.
Ldone:
z = y * (LOG2E - 1.0);
z += x * (LOG2E - 1.0);
z += y;
z += x;
z += exp;
return z;
}
}
///
unittest
{
// check if values are equal to 19 decimal digits of precision
assert(equalsDigit(log2(1024.0L), 10, 19));
}
/*****************************************
* Extracts the exponent of x as a signed integral value.
*
* If x is subnormal, it is treated as if it were normalized.
* For a positive, finite x:
*
* 1 $(LT)= $(I x) * FLT_RADIX$(SUPERSCRIPT -logb(x)) $(LT) FLT_RADIX
*
* $(TABLE_SV
* $(TR $(TH x) $(TH logb(x)) $(TH divide by 0?) )
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
* $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) )
* )
*/
real logb(real x) @trusted nothrow @nogc
{
version (Win64_DMD_InlineAsm)
{
asm pure nothrow @nogc
{
naked ;
fld real ptr [RCX] ;
fxtract ;
fstp ST(0) ;
ret ;
}
}
else version (CRuntime_Microsoft)
{
asm pure nothrow @nogc
{
fld x ;
fxtract ;
fstp ST(0) ;
}
}
else
return core.stdc.math.logbl(x);
}
/************************************
* Calculates the remainder from the calculation x/y.
* Returns:
* The value of x - i * y, where i is the number of times that y can
* be completely subtracted from x. The result has the same sign as x.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH y) $(TH fmod(x, y)) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD not 0.0) $(TD $(PLUSMN)0.0) $(TD no))
* $(TR $(TD $(PLUSMNINF)) $(TD anything) $(TD $(NAN)) $(TD yes))
* $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD $(NAN)) $(TD yes))
* $(TR $(TD !=$(PLUSMNINF)) $(TD $(PLUSMNINF)) $(TD x) $(TD no))
* )
*/
real fmod(real x, real y) @trusted nothrow @nogc
{
version (CRuntime_Microsoft)
{
return x % y;
}
else
return core.stdc.math.fmodl(x, y);
}
/************************************
* Breaks x into an integral part and a fractional part, each of which has
* the same sign as x. The integral part is stored in i.
* Returns:
* The fractional part of x.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH i (on input)) $(TH modf(x, i)) $(TH i (on return)))
* $(TR $(TD $(PLUSMNINF)) $(TD anything) $(TD $(PLUSMN)0.0) $(TD $(PLUSMNINF)))
* )
*/
real modf(real x, ref real i) @trusted nothrow @nogc
{
version (CRuntime_Microsoft)
{
i = trunc(x);
return copysign(isInfinity(x) ? 0.0 : x - i, x);
}
else
return core.stdc.math.modfl(x,&i);
}
/*************************************
* Efficiently calculates x * 2$(SUPERSCRIPT n).
*
* scalbn handles underflow and overflow in
* the same fashion as the basic arithmetic operators.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH scalb(x)))
* $(TR $(TD $(PLUSMNINF)) $(TD $(PLUSMNINF)) )
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
* )
*/
real scalbn(real x, int n) @trusted nothrow @nogc
{
version(InlineAsm_X86_Any) {
// scalbnl is not supported on DMD-Windows, so use asm pure nothrow @nogc.
version (Win64)
{
asm pure nothrow @nogc {
naked ;
mov 16[RSP],RCX ;
fild word ptr 16[RSP] ;
fld real ptr [RDX] ;
fscale ;
fstp ST(1) ;
ret ;
}
}
else
{
asm pure nothrow @nogc {
fild n;
fld x;
fscale;
fstp ST(1);
}
}
}
else
{
return core.stdc.math.scalbnl(x, n);
}
}
///
@safe nothrow @nogc unittest
{
assert(scalbn(-real.infinity, 5) == -real.infinity);
}
/***************
* Calculates the cube root of x.
*
* $(TABLE_SV
* $(TR $(TH $(I x)) $(TH cbrt(x)) $(TH invalid?))
* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) )
* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) )
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) )
* )
*/
real cbrt(real x) @trusted nothrow @nogc
{
version (CRuntime_Microsoft)
{
version (INLINE_YL2X)
return copysign(exp2(yl2x(fabs(x), 1.0L/3.0L)), x);
else
return core.stdc.math.cbrtl(x);
}
else
return core.stdc.math.cbrtl(x);
}
/*******************************
* Returns |x|
*
* $(TABLE_SV
* $(TR $(TH x) $(TH fabs(x)))
* $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) )
* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) )
* )
*/
real fabs(real x) @safe pure nothrow @nogc; /* intrinsic */
//FIXME
///ditto
double fabs(double x) @safe pure nothrow @nogc { return fabs(cast(real)x); }
//FIXME
///ditto
float fabs(float x) @safe pure nothrow @nogc { return fabs(cast(real)x); }
/***********************************************************************
* Calculates the length of the
* hypotenuse of a right-angled triangle with sides of length x and y.
* The hypotenuse is the value of the square root of
* the sums of the squares of x and y:
*
* sqrt($(POWER x, 2) + $(POWER y, 2))
*
* Note that hypot(x, y), hypot(y, x) and
* hypot(x, -y) are equivalent.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH y) $(TH hypot(x, y)) $(TH invalid?))
* $(TR $(TD x) $(TD $(PLUSMN)0.0) $(TD |x|) $(TD no))
* $(TR $(TD $(PLUSMNINF)) $(TD y) $(TD +$(INFIN)) $(TD no))
* $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD +$(INFIN)) $(TD no))
* )
*/
real hypot(real x, real y) @safe pure nothrow @nogc
{
// Scale x and y to avoid underflow and overflow.
// If one is huge and the other tiny, return the larger.
// If both are huge, avoid overflow by scaling by 1/sqrt(real.max/2).
// If both are tiny, avoid underflow by scaling by sqrt(real.min_normal*real.epsilon).
enum real SQRTMIN = 0.5 * sqrt(real.min_normal); // This is a power of 2.
enum real SQRTMAX = 1.0L / SQRTMIN; // 2^^((max_exp)/2) = nextUp(sqrt(real.max))
static assert(2*(SQRTMAX/2)*(SQRTMAX/2) <= real.max);
// Proves that sqrt(real.max) ~~ 0.5/sqrt(real.min_normal)
static assert(real.min_normal*real.max > 2 && real.min_normal*real.max <= 4);
real u = fabs(x);
real v = fabs(y);
if (!(u >= v)) // check for NaN as well.
{
v = u;
u = fabs(y);
if (u == real.infinity) return u; // hypot(inf, nan) == inf
if (v == real.infinity) return v; // hypot(nan, inf) == inf
}
// Now u >= v, or else one is NaN.
if (v >= SQRTMAX*0.5)
{
// hypot(huge, huge) -- avoid overflow
u *= SQRTMIN*0.5;
v *= SQRTMIN*0.5;
return sqrt(u*u + v*v) * SQRTMAX * 2.0;
}
if (u <= SQRTMIN)
{
// hypot (tiny, tiny) -- avoid underflow
// This is only necessary to avoid setting the underflow
// flag.
u *= SQRTMAX / real.epsilon;
v *= SQRTMAX / real.epsilon;
return sqrt(u*u + v*v) * SQRTMIN * real.epsilon;
}
if (u * real.epsilon > v)
{
// hypot (huge, tiny) = huge
return u;
}
// both are in the normal range
return sqrt(u*u + v*v);
}
unittest
{
static real[3][] vals = // x,y,hypot
[
[ 0.0, 0.0, 0.0],
[ 0.0, -0.0, 0.0],
[ -0.0, -0.0, 0.0],
[ 3.0, 4.0, 5.0],
[ -300, -400, 500],
[0.0, 7.0, 7.0],
[9.0, 9*real.epsilon, 9.0],
[88/(64*sqrt(real.min_normal)), 105/(64*sqrt(real.min_normal)), 137/(64*sqrt(real.min_normal))],
[88/(128*sqrt(real.min_normal)), 105/(128*sqrt(real.min_normal)), 137/(128*sqrt(real.min_normal))],
[3*real.min_normal*real.epsilon, 4*real.min_normal*real.epsilon, 5*real.min_normal*real.epsilon],
[ real.min_normal, real.min_normal, sqrt(2.0L)*real.min_normal],
[ real.max/sqrt(2.0L), real.max/sqrt(2.0L), real.max],
[ real.infinity, real.nan, real.infinity],
[ real.nan, real.infinity, real.infinity],
[ real.nan, real.nan, real.nan],
[ real.nan, real.max, real.nan],
[ real.max, real.nan, real.nan],
];
for (int i = 0; i < vals.length; i++)
{
real x = vals[i][0];
real y = vals[i][1];
real z = vals[i][2];
real h = hypot(x, y);
assert(isIdentical(z,h) || feqrel(z, h) >= real.mant_dig - 1);
}
}
/**************************************
* Returns the value of x rounded upward to the next integer
* (toward positive infinity).
*/
real ceil(real x) @trusted pure nothrow @nogc
{
version (Win64_DMD_InlineAsm)
{
asm pure nothrow @nogc
{
naked ;
fld real ptr [RCX] ;
fstcw 8[RSP] ;
mov AL,9[RSP] ;
mov DL,AL ;
and AL,0xC3 ;
or AL,0x08 ; // round to +infinity
mov 9[RSP],AL ;
fldcw 8[RSP] ;
frndint ;
mov 9[RSP],DL ;
fldcw 8[RSP] ;
ret ;
}
}
else version(CRuntime_Microsoft)
{
short cw;
asm pure nothrow @nogc
{
fld x ;
fstcw cw ;
mov AL,byte ptr cw+1 ;
mov DL,AL ;
and AL,0xC3 ;
or AL,0x08 ; // round to +infinity
mov byte ptr cw+1,AL ;
fldcw cw ;
frndint ;
mov byte ptr cw+1,DL ;
fldcw cw ;
}
}
else
{
// Special cases.
if (isNaN(x) || isInfinity(x))
return x;
real y = floorImpl(x);
if (y < x)
y += 1.0;
return y;
}
}
///
@safe pure nothrow @nogc unittest
{
assert(ceil(+123.456L) == +124);
assert(ceil(-123.456L) == -123);
assert(ceil(-1.234L) == -1);
assert(ceil(-0.123L) == 0);
assert(ceil(0.0L) == 0);
assert(ceil(+0.123L) == 1);
assert(ceil(+1.234L) == 2);
assert(ceil(real.infinity) == real.infinity);
assert(isNaN(ceil(real.nan)));
assert(isNaN(ceil(real.init)));
}
// ditto
double ceil(double x) @trusted pure nothrow @nogc
{
// Special cases.
if (isNaN(x) || isInfinity(x))
return x;
double y = floorImpl(x);
if (y < x)
y += 1.0;
return y;
}
@safe pure nothrow @nogc unittest
{
assert(ceil(+123.456) == +124);
assert(ceil(-123.456) == -123);
assert(ceil(-1.234) == -1);
assert(ceil(-0.123) == 0);
assert(ceil(0.0) == 0);
assert(ceil(+0.123) == 1);
assert(ceil(+1.234) == 2);
assert(ceil(double.infinity) == double.infinity);
assert(isNaN(ceil(double.nan)));
assert(isNaN(ceil(double.init)));
}
// ditto
float ceil(float x) @trusted pure nothrow @nogc
{
// Special cases.
if (isNaN(x) || isInfinity(x))
return x;
float y = floorImpl(x);
if (y < x)
y += 1.0;
return y;
}
@safe pure nothrow @nogc unittest
{
assert(ceil(+123.456f) == +124);
assert(ceil(-123.456f) == -123);
assert(ceil(-1.234f) == -1);
assert(ceil(-0.123f) == 0);
assert(ceil(0.0f) == 0);
assert(ceil(+0.123f) == 1);
assert(ceil(+1.234f) == 2);
assert(ceil(float.infinity) == float.infinity);
assert(isNaN(ceil(float.nan)));
assert(isNaN(ceil(float.init)));
}
/**************************************
* Returns the value of x rounded downward to the next integer
* (toward negative infinity).
*/
real floor(real x) @trusted pure nothrow @nogc
{
version (Win64_DMD_InlineAsm)
{
asm pure nothrow @nogc
{
naked ;
fld real ptr [RCX] ;
fstcw 8[RSP] ;
mov AL,9[RSP] ;
mov DL,AL ;
and AL,0xC3 ;
or AL,0x04 ; // round to -infinity
mov 9[RSP],AL ;
fldcw 8[RSP] ;
frndint