aboutsummaryrefslogtreecommitdiff
path: root/libphobos/src/std/math.d
diff options
context:
space:
mode:
Diffstat (limited to 'libphobos/src/std/math.d')
-rw-r--r--libphobos/src/std/math.d8413
1 files changed, 8413 insertions, 0 deletions
diff --git a/libphobos/src/std/math.d b/libphobos/src/std/math.d
new file mode 100644
index 0000000..84cf4a7
--- /dev/null
+++ b/libphobos/src/std/math.d
@@ -0,0 +1,8413 @@
+// 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) $(MYREF nextPow2) $(MYREF truncPow2)
+))
+$(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) $(MYREF quantize)
+))
+$(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)
+ $(MYREF cmp)
+))
+$(TR $(TDNW Introspection) $(TD
+ $(MYREF isFinite) $(MYREF isIdentical) $(MYREF isInfinity) $(MYREF isNaN)
+ $(MYREF isNormal) $(MYREF isSubnormal) $(MYREF signbit) $(MYREF sgn)
+ $(MYREF copysign) $(MYREF isPowerOf2)
+))
+$(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:
+ * 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))
+ * TABLE_DOMRG = <table border="1" cellpadding="4" cellspacing="0">
+ * $(SVH Domain X, Range Y)
+ $(SV $1, $2)
+ * </table>
+ * DOMAIN=$1
+ * RANGE=$1
+
+ * 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 $(LT)steve@moshier.net$(GT)
+ * 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: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
+ * Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston,
+ * Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
+ * Source: $(PHOBOSSRC std/_math.d)
+ */
+
+/* NOTE: This file has been patched from the original DMD distribution to
+ * work with the GDC compiler.
+ */
+module std.math;
+
+version (Win64)
+{
+ version (D_InlineAsm_X86_64)
+ version = Win64_DMD_InlineAsm;
+}
+
+static import core.math;
+static import core.stdc.math;
+static import core.stdc.fenv;
+import std.traits; // CommonType, isFloatingPoint, isIntegral, isSigned, isUnsigned, Largest, Unqual
+
+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 (MIPS32) version = MIPS_Any;
+version (MIPS64) version = MIPS_Any;
+version (AArch64) version = ARM_Any;
+version (ARM) version = ARM_Any;
+
+version (D_InlineAsm_X86)
+{
+ version = InlineAsm_X86_Any;
+}
+else version (D_InlineAsm_X86_64)
+{
+ version = InlineAsm_X86_Any;
+}
+
+version (X86_64) version = StaticallyHaveSSE;
+version (X86) version (OSX) version = StaticallyHaveSSE;
+
+version (StaticallyHaveSSE)
+{
+ private enum bool haveSSE = true;
+}
+else
+{
+ static import core.cpuid;
+ private alias haveSSE = core.cpuid.sse;
+}
+
+version (unittest)
+{
+ import core.stdc.stdio; // : sprintf;
+
+ 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];
+ }
+}
+
+
+
+package:
+// 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
+ // EXPBIAS is the exponent bias - 1 (exp == EXPBIAS yields ×2^-1).
+ // 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 = 0x3FFE;
+ 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.meta : AliasSeq;
+ foreach (T; AliasSeq!(float, double, real))
+ {
+ T f = 3;
+ assert(abs(f) == f);
+ assert(abs(-f) == f);
+ }
+ foreach (T; AliasSeq!(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 { pragma(inline, true); return core.math.cos(x); }
+//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); }
+
+@safe unittest
+{
+ real function(real) pcos = &cos;
+ assert(pcos != null);
+}
+
+/***********************************
+ * Returns $(HTTP en.wikipedia.org/wiki/Sine, sine) of x. x is in $(HTTP 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 { pragma(inline, true); return core.math.sin(x); }
+//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); }
+
+///
+@safe 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);
+ }
+}
+
+@safe unittest
+{
+ real function(real) psin = &sin;
+ assert(psin != null);
+}
+
+/***********************************
+ * 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
+{
+ const creal cs = expi(z.re);
+ const 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
+{
+ const creal cs = expi(z.re);
+ const 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 ;
+ fstsw AX ;
+ sahf ;
+ jnp Clear1 ; // 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;
+
+Clear1: asm pure nothrow @nogc{
+ fstp ST(0) ; // dump X, which is always 1
+ }
+
+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 ;
+ fstsw AX ;
+ test AH,4 ;
+ jz Clear1 ; // 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;
+
+Clear1: asm pure nothrow @nogc{
+ fstp ST(0) ; // dump X, which is always 1
+ }
+
+Lret: {}
+ }
+ else
+ {
+ // Coefficients for tan(x) and PI/4 split into three parts.
+ static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple)
+ {
+ static immutable real[6] P = [
+ 2.883414728874239697964612246732416606301E10L,
+ -2.307030822693734879744223131873392503321E9L,
+ 5.160188250214037865511600561074819366815E7L,
+ -4.249691853501233575668486667664718192660E5L,
+ 1.272297782199996882828849455156962260810E3L,
+ -9.889929415807650724957118893791829849557E-1L
+ ];
+ static immutable real[7] Q = [
+ 8.650244186622719093893836740197250197602E10L
+ -4.152206921457208101480801635640958361612E10L,
+ 2.758476078803232151774723646710890525496E9L,
+ -5.733709132766856723608447733926138506824E7L,
+ 4.529422062441341616231663543669583527923E5L,
+ -1.317243702830553658702531997959756728291E3L,
+ 1.0
+ ];
+
+ enum real P1 =
+ 7.853981633974483067550664827649598009884357452392578125E-1L;
+ enum real P2 =
+ 2.8605943630549158983813312792950660807511260829685741796657E-18L;
+ enum real P3 =
+ 2.1679525325309452561992610065108379921905808E-35L;
+ }
+ else
+ {
+ 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,
+ ];
+
+ 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;
+ const 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) ));
+}
+
+@system 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); }
+
+@system 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); }
+
+@system 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 if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple)
+ {
+ static immutable real[9] P = [
+ -6.880597774405940432145577545328795037141E2L,
+ -2.514829758941713674909996882101723647996E3L,
+ -3.696264445691821235400930243493001671932E3L,
+ -2.792272753241044941703278827346430350236E3L,
+ -1.148164399808514330375280133523543970854E3L,
+ -2.497759878476618348858065206895055957104E2L,
+ -2.548067867495502632615671450650071218995E1L,
+ -8.768423468036849091777415076702113400070E-1L,
+ -6.635810778635296712545011270011752799963E-4L
+ ];
+ static immutable real[9] Q = [
+ 2.064179332321782129643673263598686441900E3L,
+ 8.782996876218210302516194604424986107121E3L,
+ 1.547394317752562611786521896296215170819E4L,
+ 1.458510242529987155225086911411015961174E4L,
+ 7.928572347062145288093560392463784743935E3L,
+ 2.494680540950601626662048893678584497900E3L,
+ 4.308348370818927353321556740027020068897E2L,
+ 3.566239794444800849656497338030115886153E1L,
+ 1.0
+ ];
+ }
+ else
+ {
+ 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 = 0.414213562373095048801688724209698078569672L;
+ // tan(3 * PI/8)
+ enum real TAN3_PI_8 = 2.414213562373095048801688724209698078569672L;
+
+ // 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.
+ const 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); }
+
+@system 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);
+}
+
+@system 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.
+ const 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); }
+
+@system 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);
+ }
+
+ 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); }
+
+@system 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);
+ }
+
+ const 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); }
+
+@system 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)
+ {
+ const real y = exp(fabs(x));
+ return y * 0.5 + 0.5i * copysign(y, x);
+ }
+ else
+ {
+ const 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 0..$(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); }
+
+
+@system 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); }
+
+@system 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)
+ * )
+ * $(BR)
+ * $(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); }
+
+
+@system 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 { pragma(inline, true); return core.math.rndtol(x); }
+//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); }
+
+@safe unittest
+{
+ long function(real) prndtol = &rndtol;
+ assert(prndtol != null);
+}
+
+/*****************************************
+ * 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 { pragma(inline, true); return core.math.sqrt(x); }
+
+/// ditto
+double sqrt(double x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.sqrt(x); }
+
+/// ditto
+real sqrt(real x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.sqrt(x); }
+
+@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)));
+}
+
+@safe unittest
+{
+ float function(float) psqrtf = &sqrt;
+ assert(psqrtf != null);
+ double function(double) psqrtd = &sqrt;
+ assert(psqrtd != null);
+ real function(real) psqrtr = &sqrt;
+ assert(psqrtr != null);
+}
+
+creal sqrt(creal z) @nogc @safe pure nothrow
+{
+ creal c;
+ real x,y,w,r;
+
+ if (z == 0)
+ {
+ c = 0 + 0i;
+ }
+ else
+ {
+ const real z_re = z.re;
+ const 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
+ {
+ alias F = floatTraits!real;
+ static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ // Coefficients for exp(x)
+ static immutable real[3] P = [
+ 9.99999999999999999910E-1L,
+ 3.02994407707441961300E-2L,
+ 1.26177193074810590878E-4L,
+ ];
+ static immutable real[4] Q = [
+ 2.00000000000000000009E0L,
+ 2.27265548208155028766E-1L,
+ 2.52448340349684104192E-3L,
+ 3.00198505138664455042E-6L,
+ ];
+
+ // C1 + C2 = LN2.
+ enum real C1 = 6.93145751953125E-1;
+ enum real C2 = 1.42860682030941723212E-6;
+
+ // Overflow and Underflow limits.
+ enum real OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024)
+ enum real UF = -7.451332191019412076235E2; // ln(2^-1075)
+ }
+ else static if (F.realFormat == RealFormat.ieeeExtended)
+ {
+ // 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.4286068203094172321215E-6L;
+
+ // Overflow and Underflow limits.
+ enum real OF = 1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384)
+ enum real UF = -1.13994985314888605586758E4L; // ln(2^-16446)
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ // Coefficients for exp(x) - 1
+ static immutable real[5] P = [
+ 9.999999999999999999999999999999999998502E-1L,
+ 3.508710990737834361215404761139478627390E-2L,
+ 2.708775201978218837374512615596512792224E-4L,
+ 6.141506007208645008909088812338454698548E-7L,
+ 3.279723985560247033712687707263393506266E-10L
+ ];
+ static immutable real[6] Q = [
+ 2.000000000000000000000000000000000000150E0,
+ 2.368408864814233538909747618894558968880E-1L,
+ 3.611828913847589925056132680618007270344E-3L,
+ 1.504792651814944826817779302637284053660E-5L,
+ 1.771372078166251484503904874657985291164E-8L,
+ 2.980756652081995192255342779918052538681E-12L
+ ];
+
+ // C1 + C2 = LN2.
+ enum real C1 = 6.93145751953125E-1L;
+ enum real C2 = 1.428606820309417232121458176568075500134E-6L;
+
+ // Overflow and Underflow limits.
+ enum real OF = 1.135583025911358400418251384584930671458833e4L;
+ enum real UF = -1.143276959615573793352782661133116431383730e4L;
+ }
+ else
+ static assert(0, "Not implemented for this architecture");
+
+ // Special cases. Raises an overflow or underflow flag accordingly,
+ // except in the case for CTFE, where there are no hardware controls.
+ if (isNaN(x))
+ return x;
+ if (x > OF)
+ {
+ if (__ctfe)
+ return real.infinity;
+ else
+ return real.max * copysign(real.max, real.infinity);
+ }
+ if (x < UF)
+ {
+ if (__ctfe)
+ return 0.0;
+ else
+ return real.min_normal * copysign(real.min_normal, 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))
+ const real xx = x * x;
+ const 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); }
+
+@system 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 and overflow/underflow limits.
+ static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple)
+ {
+ static immutable real[8] P = [
+ 2.943520915569954073888921213330863757240E8L,
+ -5.722847283900608941516165725053359168840E7L,
+ 8.944630806357575461578107295909719817253E6L,
+ -7.212432713558031519943281748462837065308E5L,
+ 4.578962475841642634225390068461943438441E4L,
+ -1.716772506388927649032068540558788106762E3L,
+ 4.401308817383362136048032038528753151144E1L,
+ -4.888737542888633647784737721812546636240E-1L
+ ];
+
+ static immutable real[9] Q = [
+ 1.766112549341972444333352727998584753865E9L,
+ -7.848989743695296475743081255027098295771E8L,
+ 1.615869009634292424463780387327037251069E8L,
+ -2.019684072836541751428967854947019415698E7L,
+ 1.682912729190313538934190635536631941751E6L,
+ -9.615511549171441430850103489315371768998E4L,
+ 3.697714952261803935521187272204485251835E3L,
+ -8.802340681794263968892934703309274564037E1L,
+ 1.0
+ ];
+
+ enum real OF = 1.1356523406294143949491931077970764891253E4L;
+ enum real UF = -1.143276959615573793352782661133116431383730e4L;
+ }
+ else
+ {
+ 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.0
+ ];
+
+ enum real OF = 1.1356523406294143949492E4L;
+ enum real UF = -4.5054566736396445112120088E1L;
+ }
+
+
+ // C1 + C2 = LN2.
+ enum real C1 = 6.9314575195312500000000E-1L;
+ enum real C2 = 1.428606820309417232121458176568075500134E-6L;
+
+ // Special cases. Raises an overflow flag, except in the case
+ // for CTFE, where there are no hardware controls.
+ if (x > OF)
+ {
+ if (__ctfe)
+ return real.infinity;
+ else
+ return real.max * copysign(real.max, 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);
+ const 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)) )
+ * )
+ */
+pragma(inline, true)
+real exp2(real x) @nogc @trusted pure nothrow
+{
+ version (InlineAsm_X86_Any)
+ {
+ if (!__ctfe)
+ return exp2Asm(x);
+ else
+ return exp2Impl(x);
+ }
+ else
+ {
+ return exp2Impl(x);
+ }
+}
+
+version (InlineAsm_X86_Any)
+private real exp2Asm(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
+ static assert(0);
+}
+
+private real exp2Impl(real x) @nogc @trusted pure nothrow
+{
+ // Coefficients for exp2(x)
+ static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple)
+ {
+ static immutable real[5] P = [
+ 9.079594442980146270952372234833529694788E12L,
+ 1.530625323728429161131811299626419117557E11L,
+ 5.677513871931844661829755443994214173883E8L,
+ 6.185032670011643762127954396427045467506E5L,
+ 1.587171580015525194694938306936721666031E2L
+ ];
+
+ static immutable real[6] Q = [
+ 2.619817175234089411411070339065679229869E13L,
+ 1.490560994263653042761789432690793026977E12L,
+ 1.092141473886177435056423606755843616331E10L,
+ 2.186249607051644894762167991800811827835E7L,
+ 1.236602014442099053716561665053645270207E4L,
+ 1.0
+ ];
+ }
+ else
+ {
+ 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 = 16_384.0L;
+ enum real UF = -16_382.0L;
+
+ // Special cases. Raises an overflow or underflow flag accordingly,
+ // except in the case for CTFE, where there are no hardware controls.
+ if (isNaN(x))
+ return x;
+ if (x > OF)
+ {
+ if (__ctfe)
+ return real.infinity;
+ else
+ return real.max * copysign(real.max, real.infinity);
+ }
+ if (x < UF)
+ {
+ if (__ctfe)
+ return 0.0;
+ else
+ return real.min_normal * copysign(real.min_normal, 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))
+ const real xx = x * x;
+ const 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;
+}
+
+///
+@safe 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);
+}
+
+@safe 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 );
+ }
+}
+
+@system unittest
+{
+ FloatingPointControl ctrl;
+ if (FloatingPointControl.hasExceptionTraps)
+ ctrl.disableExceptions(FloatingPointControl.allExceptions);
+ ctrl.rounding = FloatingPointControl.roundToNearest;
+
+ static if (real.mant_dig == 113)
+ {
+ static immutable real[2][] exptestpoints =
+ [ // x exp(x)
+ [ 1.0L, E ],
+ [ 0.5L, 0x1.a61298e1e069bc972dfefab6df34p+0L ],
+ [ 3.0L, E*E*E ],
+ [ 0x1.6p+13L, 0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow
+ [ 0x1.7p+13L, real.infinity ], // close overflow
+ [ 0x1p+80L, real.infinity ], // far overflow
+ [ real.infinity, real.infinity ],
+ [-0x1.18p+13L, 0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow
+ [-0x1.625p+13L, 0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto
+ [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal
+ [-0x1.6549p+13L, 0x0.0000000000000000000000000001p-16382L ], // ditto
+ [-0x1.655p+13L, 0 ], // close underflow
+ [-0x1p+30L, 0 ], // far underflow
+ ];
+ }
+ else 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 minEqualDecimalDigits = real.dig - 3;
+ real x;
+ IeeeFlags f;
+ foreach (ref pair; exptestpoints)
+ {
+ resetIeeeFlags();
+ x = exp(pair[0]);
+ f = ieeeFlags;
+ assert(equalsDigit(x, pair[1], minEqualDecimalDigits));
+
+ 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 (verified against GNU MPFR/Mathematica).
+ assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+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] = ((-1 - F.EXPMASK) & 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] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
+ }
+ }
+ 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] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
+ }
+ 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)(((-1 - F.EXPMASK) & 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)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00);
+ }
+ return vf;
+ }
+ else // static if (F.realFormat == RealFormat.ibmExtended)
+ {
+ assert(0, "frexp not implemented");
+ }
+}
+
+///
+@system 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);
+}
+
+@safe unittest
+{
+ import std.meta : AliasSeq;
+ import std.typecons : tuple, Tuple;
+
+ foreach (T; AliasSeq!(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),
+
+ // Phobos issue #16026:
+ tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
+ ];
+
+ 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);
+
+ }
+ }
+ }
+}
+
+@safe unittest
+{
+ import std.meta : AliasSeq;
+ void foo() {
+ foreach (T; AliasSeq!(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(T)(const T x) @trusted pure nothrow @nogc
+if (isFloatingPoint!T)
+{
+ import core.bitop : bsr;
+ alias F = floatTraits!T;
+
+ union floatBits
+ {
+ T rv;
+ ushort[T.sizeof/2] vu;
+ uint[T.sizeof/4] vui;
+ static if (T.sizeof >= 8)
+ ulong[T.sizeof/8] vul;
+ }
+ floatBits y = void;
+ y.rv = x;
+
+ int ex = y.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 (y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) // NaN
+ return FP_ILOGBNAN;
+ else // +-infinity
+ return int.max;
+ }
+ else
+ {
+ return ex - F.EXPBIAS - 1;
+ }
+ }
+ else if (!y.vul[0])
+ {
+ // vf is +-0.0
+ return FP_ILOGB0;
+ }
+ else
+ {
+ // subnormal
+ return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]);
+ }
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ if (ex) // If exponent is non-zero
+ {
+ if (ex == F.EXPMASK)
+ {
+ // infinity or NaN
+ if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN
+ return FP_ILOGBNAN;
+ else // +- infinity
+ return int.max;
+ }
+ else
+ {
+ return ex - F.EXPBIAS - 1;
+ }
+ }
+ else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
+ {
+ // vf is +-0.0
+ return FP_ILOGB0;
+ }
+ else
+ {
+ // subnormal
+ const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF;
+ const ulong lsb = y.vul[MANTISSA_LSB];
+ if (msb)
+ return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64;
+ else
+ return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb);
+ }
+ }
+ else static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ if (ex) // If exponent is non-zero
+ {
+ if (ex == F.EXPMASK) // infinity or NaN
+ {
+ if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000) // +- infinity
+ return int.max;
+ else // NaN
+ return FP_ILOGBNAN;
+ }
+ else
+ {
+ return ((ex - F.EXPBIAS) >> 4) - 1;
+ }
+ }
+ else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF))
+ {
+ // vf is +-0.0
+ return FP_ILOGB0;
+ }
+ else
+ {
+ // subnormal
+ enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF;
+ return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64);
+ }
+ }
+ else static if (F.realFormat == RealFormat.ieeeSingle)
+ {
+ if (ex) // If exponent is non-zero
+ {
+ if (ex == F.EXPMASK) // infinity or NaN
+ {
+ if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000) // +- infinity
+ return int.max;
+ else // NaN
+ return FP_ILOGBNAN;
+ }
+ else
+ {
+ return ((ex - F.EXPBIAS) >> 7) - 1;
+ }
+ }
+ else if (!(y.vui[0] & 0x7FFF_FFFF))
+ {
+ // vf is +-0.0
+ return FP_ILOGB0;
+ }
+ else
+ {
+ // subnormal
+ const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT;
+ return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa);
+ }
+ }
+ else // static if (F.realFormat == RealFormat.ibmExtended)
+ {
+ core.stdc.math.ilogbl(x);
+ }
+}
+/// ditto
+int ilogb(T)(const T x) @safe pure nothrow @nogc
+if (isIntegral!T && isUnsigned!T)
+{
+ import core.bitop : bsr;
+ if (x == 0)
+ return FP_ILOGB0;
+ else
+ {
+ static assert(T.sizeof <= ulong.sizeof, "integer size too large for the current ilogb implementation");
+ return bsr(x);
+ }
+}
+/// ditto
+int ilogb(T)(const T x) @safe pure nothrow @nogc
+if (isIntegral!T && isSigned!T)
+{
+ import std.traits : Unsigned;
+ // Note: abs(x) can not be used because the return type is not Unsigned and
+ // the return value would be wrong for x == int.min
+ Unsigned!T absx = x >= 0 ? x : -x;
+ return ilogb(absx);
+}
+
+alias FP_ILOGB0 = core.stdc.math.FP_ILOGB0;
+alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN;
+
+@system nothrow @nogc unittest
+{
+ import std.meta : AliasSeq;
+ import std.typecons : Tuple;
+ foreach (F; AliasSeq!(float, double, real))
+ {
+ alias T = Tuple!(F, int);
+ T[13] vals = // x, ilogb(x)
+ [
+ T( F.nan , FP_ILOGBNAN ),
+ T( -F.nan , FP_ILOGBNAN ),
+ T( F.infinity, int.max ),
+ T( -F.infinity, int.max ),
+ T( 0.0 , FP_ILOGB0 ),
+ T( -0.0 , FP_ILOGB0 ),
+ T( 2.0 , 1 ),
+ T( 2.0001 , 1 ),
+ T( 1.9999 , 0 ),
+ T( 0.5 , -1 ),
+ T( 123.123 , 6 ),
+ T( -123.123 , 6 ),
+ T( 0.123 , -4 ),
+ ];
+
+ foreach (elem; vals)
+ {
+ assert(ilogb(elem[0]) == elem[1]);
+ }
+ }
+
+ // min_normal and subnormals
+ assert(ilogb(-float.min_normal) == -126);
+ assert(ilogb(nextUp(-float.min_normal)) == -127);
+ assert(ilogb(nextUp(-float(0.0))) == -149);
+ assert(ilogb(-double.min_normal) == -1022);
+ assert(ilogb(nextUp(-double.min_normal)) == -1023);
+ assert(ilogb(nextUp(-double(0.0))) == -1074);
+ static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
+ {
+ assert(ilogb(-real.min_normal) == -16382);
+ assert(ilogb(nextUp(-real.min_normal)) == -16383);
+ assert(ilogb(nextUp(-real(0.0))) == -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(-real(0.0))) == -1074);
+ }
+
+ // test integer types
+ assert(ilogb(0) == FP_ILOGB0);
+ assert(ilogb(int.max) == 30);
+ assert(ilogb(int.min) == 31);
+ assert(ilogb(uint.max) == 31);
+ assert(ilogb(long.max) == 62);
+ assert(ilogb(long.min) == 63);
+ assert(ilogb(ulong.max) == 63);
+}
+
+/*******************************************
+ * Compute n * 2$(SUPERSCRIPT exp)
+ * References: frexp
+ */
+
+real ldexp(real n, int exp) @nogc @safe pure nothrow { pragma(inline, true); return core.math.ldexp(n, exp); }
+//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.meta : AliasSeq;
+ foreach (T; AliasSeq!(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");
+}
+
+/* workaround Issue 14718, float parsing depends on platform strtold
+@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);
+}
+*/
+
+@system 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));
+ }
+
+ real function(real, int) pldexp = &ldexp;
+ assert(pldexp != null);
+}
+
+private
+{
+ version (INLINE_YL2X) {} else
+ {
+ static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple)
+ {
+ // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x)
+ static immutable real[13] logCoeffsP = [
+ 1.313572404063446165910279910527789794488E4L,
+ 7.771154681358524243729929227226708890930E4L,
+ 2.014652742082537582487669938141683759923E5L,
+ 3.007007295140399532324943111654767187848E5L,
+ 2.854829159639697837788887080758954924001E5L,
+ 1.797628303815655343403735250238293741397E5L,
+ 7.594356839258970405033155585486712125861E4L,
+ 2.128857716871515081352991964243375186031E4L,
+ 3.824952356185897735160588078446136783779E3L,
+ 4.114517881637811823002128927449878962058E2L,
+ 2.321125933898420063925789532045674660756E1L,
+ 4.998469661968096229986658302195402690910E-1L,
+ 1.538612243596254322971797716843006400388E-6L
+ ];
+ static immutable real[13] logCoeffsQ = [
+ 3.940717212190338497730839731583397586124E4L,
+ 2.626900195321832660448791748036714883242E5L,
+ 7.777690340007566932935753241556479363645E5L,
+ 1.347518538384329112529391120390701166528E6L,
+ 1.514882452993549494932585972882995548426E6L,
+ 1.158019977462989115839826904108208787040E6L,
+ 6.132189329546557743179177159925690841200E5L,
+ 2.248234257620569139969141618556349415120E5L,
+ 5.605842085972455027590989944010492125825E4L,
+ 9.147150349299596453976674231612674085381E3L,
+ 9.104928120962988414618126155557301584078E2L,
+ 4.839208193348159620282142911143429644326E1L,
+ 1.0
+ ];
+
+ // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2)
+ // where z = 2(x-1)/(x+1)
+ static immutable real[6] logCoeffsR = [
+ -8.828896441624934385266096344596648080902E-1L,
+ 8.057002716646055371965756206836056074715E1L,
+ -2.024301798136027039250415126250455056397E3L,
+ 2.048819892795278657810231591630928516206E4L,
+ -8.977257995689735303686582344659576526998E4L,
+ 1.418134209872192732479751274970992665513E5L
+ ];
+ static immutable real[6] logCoeffsS = [
+ 1.701761051846631278975701529965589676574E6L
+ -1.332535117259762928288745111081235577029E6L,
+ 4.001557694070773974936904547424676279307E5L,
+ -5.748542087379434595104154610899551484314E4L,
+ 3.998526750980007367835804959888064681098E3L,
+ -1.186359407982897997337150403816839480438E2L,
+ 1.0
+ ];
+ }
+ else
+ {
+ // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x)
+ static immutable real[7] logCoeffsP = [
+ 2.0039553499201281259648E1L,
+ 5.7112963590585538103336E1L,
+ 6.0949667980987787057556E1L,
+ 2.9911919328553073277375E1L,
+ 6.5787325942061044846969E0L,
+ 4.9854102823193375972212E-1L,
+ 4.5270000862445199635215E-5L,
+ ];
+ static immutable real[7] logCoeffsQ = [
+ 6.0118660497603843919306E1L,
+ 2.1642788614495947685003E2L,
+ 3.0909872225312059774938E2L,
+ 2.2176239823732856465394E2L,
+ 8.3047565967967209469434E1L,
+ 1.5062909083469192043167E1L,
+ 1.0000000000000000000000E0L,
+ ];
+
+ // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2)
+ // where z = 2(x-1)/(x+1)
+ static immutable real[4] logCoeffsR = [
+ -3.5717684488096787370998E1L,
+ 1.0777257190312272158094E1L,
+ -7.1990767473014147232598E-1L,
+ 1.9757429581415468984296E-3L,
+ ];
+ static immutable real[4] logCoeffsS = [
+ -4.2861221385716144629696E2L,
+ 1.9361891836232102174846E2L,
+ -2.6201045551331104417768E1L,
+ 1.0000000000000000000000E0L,
+ ];
+ }
+ }
+}
+
+/**************************************
+ * 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 core.math.yl2x(x, LN2);
+ else
+ {
+ // C1 + C2 = LN2.
+ enum real C1 = 6.93145751953125E-1L;
+ enum real C2 = 1.428606820309417232121458176568075500134E-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 R(z) / S(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, logCoeffsR) / poly(z, logCoeffsS));
+ 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, logCoeffsP) / poly(x, logCoeffsQ));
+ 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 core.math.yl2x(x, LOG2);
+ else
+ {
+ // 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 R(z) / S(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, logCoeffsR) / poly(z, logCoeffsS));
+ 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, logCoeffsP) / poly(x, logCoeffsQ));
+ 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) ? core.math.yl2xp1(x, LN2) : core.math.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 core.math.yl2x(x, 1);
+ else
+ {
+ // 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 R(z) / S(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, logCoeffsR) / poly(z, logCoeffsS));
+ 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, logCoeffsP) / poly(x, logCoeffsQ));
+ 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;
+ }
+}
+
+///
+@system 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(core.math.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 { pragma(inline, true); return core.math.fabs(x); }
+//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); }
+
+@safe unittest
+{
+ real function(real) pfabs = &fabs;
+ assert(pfabs != null);
+}
+
+/***********************************************************************
+ * 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);
+}
+
+@safe 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 ;
+ 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,0x04 ; // 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) || x == 0.0)
+ return x;
+
+ return floorImpl(x);
+ }
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ assert(floor(+123.456L) == +123);
+ assert(floor(-123.456L) == -124);
+ assert(floor(-1.234L) == -2);
+ assert(floor(-0.123L) == -1);
+ assert(floor(0.0L) == 0);
+ assert(floor(+0.123L) == 0);
+ assert(floor(+1.234L) == 1);
+ assert(floor(real.infinity) == real.infinity);
+ assert(isNaN(floor(real.nan)));
+ assert(isNaN(floor(real.init)));
+}
+
+// ditto
+double floor(double x) @trusted pure nothrow @nogc
+{
+ // Special cases.
+ if (isNaN(x) || isInfinity(x) || x == 0.0)
+ return x;
+
+ return floorImpl(x);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ assert(floor(+123.456) == +123);
+ assert(floor(-123.456) == -124);
+ assert(floor(-1.234) == -2);
+ assert(floor(-0.123) == -1);
+ assert(floor(0.0) == 0);
+ assert(floor(+0.123) == 0);
+ assert(floor(+1.234) == 1);
+ assert(floor(double.infinity) == double.infinity);
+ assert(isNaN(floor(double.nan)));
+ assert(isNaN(floor(double.init)));
+}
+
+// ditto
+float floor(float x) @trusted pure nothrow @nogc
+{
+ // Special cases.
+ if (isNaN(x) || isInfinity(x) || x == 0.0)
+ return x;
+
+ return floorImpl(x);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ assert(floor(+123.456f) == +123);
+ assert(floor(-123.456f) == -124);
+ assert(floor(-1.234f) == -2);
+ assert(floor(-0.123f) == -1);
+ assert(floor(0.0f) == 0);
+ assert(floor(+0.123f) == 0);
+ assert(floor(+1.234f) == 1);
+ assert(floor(float.infinity) == float.infinity);
+ assert(isNaN(floor(float.nan)));
+ assert(isNaN(floor(float.init)));
+}
+
+/**
+ * Round `val` to a multiple of `unit`. `rfunc` specifies the rounding
+ * function to use; by default this is `rint`, which uses the current
+ * rounding mode.
+ */
+Unqual!F quantize(alias rfunc = rint, F)(const F val, const F unit)
+if (is(typeof(rfunc(F.init)) : F) && isFloatingPoint!F)
+{
+ typeof(return) ret = val;
+ if (unit != 0)
+ {
+ const scaled = val / unit;
+ if (!scaled.isInfinity)
+ ret = rfunc(scaled) * unit;
+ }
+ return ret;
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ assert(12345.6789L.quantize(0.01L) == 12345.68L);
+ assert(12345.6789L.quantize!floor(0.01L) == 12345.67L);
+ assert(12345.6789L.quantize(22.0L) == 12342.0L);
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ assert(12345.6789L.quantize(0) == 12345.6789L);
+ assert(12345.6789L.quantize(real.infinity).isNaN);
+ assert(12345.6789L.quantize(real.nan).isNaN);
+ assert(real.infinity.quantize(0.01L) == real.infinity);
+ assert(real.infinity.quantize(real.nan).isNaN);
+ assert(real.nan.quantize(0.01L).isNaN);
+ assert(real.nan.quantize(real.infinity).isNaN);
+ assert(real.nan.quantize(real.nan).isNaN);
+}
+
+/**
+ * Round `val` to a multiple of `pow(base, exp)`. `rfunc` specifies the
+ * rounding function to use; by default this is `rint`, which uses the
+ * current rounding mode.
+ */
+Unqual!F quantize(real base, alias rfunc = rint, F, E)(const F val, const E exp)
+if (is(typeof(rfunc(F.init)) : F) && isFloatingPoint!F && isIntegral!E)
+{
+ // TODO: Compile-time optimization for power-of-two bases?
+ return quantize!rfunc(val, pow(cast(F) base, exp));
+}
+
+/// ditto
+Unqual!F quantize(real base, long exp = 1, alias rfunc = rint, F)(const F val)
+if (is(typeof(rfunc(F.init)) : F) && isFloatingPoint!F)
+{
+ enum unit = cast(F) pow(base, exp);
+ return quantize!rfunc(val, unit);
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ assert(12345.6789L.quantize!10(-2) == 12345.68L);
+ assert(12345.6789L.quantize!(10, -2) == 12345.68L);
+ assert(12345.6789L.quantize!(10, floor)(-2) == 12345.67L);
+ assert(12345.6789L.quantize!(10, -2, floor) == 12345.67L);
+
+ assert(12345.6789L.quantize!22(1) == 12342.0L);
+ assert(12345.6789L.quantize!22 == 12342.0L);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ import std.meta : AliasSeq;
+
+ foreach (F; AliasSeq!(real, double, float))
+ {
+ const maxL10 = cast(int) F.max.log10.floor;
+ const maxR10 = pow(cast(F) 10, maxL10);
+ assert((cast(F) 0.9L * maxR10).quantize!10(maxL10) == maxR10);
+ assert((cast(F)-0.9L * maxR10).quantize!10(maxL10) == -maxR10);
+
+ assert(F.max.quantize(F.min_normal) == F.max);
+ assert((-F.max).quantize(F.min_normal) == -F.max);
+ assert(F.min_normal.quantize(F.max) == 0);
+ assert((-F.min_normal).quantize(F.max) == 0);
+ assert(F.min_normal.quantize(F.min_normal) == F.min_normal);
+ assert((-F.min_normal).quantize(F.min_normal) == -F.min_normal);
+ }
+}
+
+/******************************************
+ * Rounds x to the nearest integer value, using the current rounding
+ * mode.
+ *
+ * Unlike the rint functions, nearbyint does not raise the
+ * FE_INEXACT exception.
+ */
+real nearbyint(real x) @trusted nothrow @nogc
+{
+ version (CRuntime_Microsoft)
+ {
+ assert(0); // not implemented in C library
+ }
+ else
+ return core.stdc.math.nearbyintl(x);
+}
+
+/**********************************
+ * Rounds x to the nearest integer value, using the current rounding
+ * mode.
+ * If the return value is not equal to x, the FE_INEXACT
+ * exception is raised.
+ * $(B nearbyint) performs
+ * the same operation, but does not set the FE_INEXACT exception.
+ */
+real rint(real x) @safe pure nothrow @nogc { pragma(inline, true); return core.math.rint(x); }
+//FIXME
+///ditto
+double rint(double x) @safe pure nothrow @nogc { return rint(cast(real) x); }
+//FIXME
+///ditto
+float rint(float x) @safe pure nothrow @nogc { return rint(cast(real) x); }
+
+@safe unittest
+{
+ real function(real) print = &rint;
+ assert(print != null);
+}
+
+/***************************************
+ * Rounds x to the nearest integer value, using the current rounding
+ * mode.
+ *
+ * This is generally the fastest method to convert a floating-point number
+ * to an integer. Note that the results from this function
+ * depend on the rounding mode, if the fractional part of x is exactly 0.5.
+ * If using the default rounding mode (ties round to even integers)
+ * lrint(4.5) == 4, lrint(5.5)==6.
+ */
+long lrint(real x) @trusted pure nothrow @nogc
+{
+ version (InlineAsm_X86_Any)
+ {
+ version (Win64)
+ {
+ asm pure nothrow @nogc
+ {
+ naked;
+ fld real ptr [RCX];
+ fistp qword ptr 8[RSP];
+ mov RAX,8[RSP];
+ ret;
+ }
+ }
+ else
+ {
+ long n;
+ asm pure nothrow @nogc
+ {
+ fld x;
+ fistp n;
+ }
+ return n;
+ }
+ }
+ else
+ {
+ alias F = floatTraits!(real);
+ static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ long result;
+
+ // Rounding limit when casting from real(double) to ulong.
+ enum real OF = 4.50359962737049600000E15L;
+
+ uint* vi = cast(uint*)(&x);
+
+ // Find the exponent and sign
+ uint msb = vi[MANTISSA_MSB];
+ uint lsb = vi[MANTISSA_LSB];
+ int exp = ((msb >> 20) & 0x7ff) - 0x3ff;
+ const int sign = msb >> 31;
+ msb &= 0xfffff;
+ msb |= 0x100000;
+
+ if (exp < 63)
+ {
+ if (exp >= 52)
+ result = (cast(long) msb << (exp - 20)) | (lsb << (exp - 52));
+ else
+ {
+ // Adjust x and check result.
+ const real j = sign ? -OF : OF;
+ x = (j + x) - j;
+ msb = vi[MANTISSA_MSB];
+ lsb = vi[MANTISSA_LSB];
+ exp = ((msb >> 20) & 0x7ff) - 0x3ff;
+ msb &= 0xfffff;
+ msb |= 0x100000;
+
+ if (exp < 0)
+ result = 0;
+ else if (exp < 20)
+ result = cast(long) msb >> (20 - exp);
+ else if (exp == 20)
+ result = cast(long) msb;
+ else
+ result = (cast(long) msb << (exp - 20)) | (lsb >> (52 - exp));
+ }
+ }
+ else
+ {
+ // It is left implementation defined when the number is too large.
+ return cast(long) x;
+ }
+
+ return sign ? -result : result;
+ }
+ else static if (F.realFormat == RealFormat.ieeeExtended)
+ {
+ long result;
+
+ // Rounding limit when casting from real(80-bit) to ulong.
+ enum real OF = 9.22337203685477580800E18L;
+
+ ushort* vu = cast(ushort*)(&x);
+ uint* vi = cast(uint*)(&x);
+
+ // Find the exponent and sign
+ int exp = (vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff;
+ const int sign = (vu[F.EXPPOS_SHORT] >> 15) & 1;
+
+ if (exp < 63)
+ {
+ // Adjust x and check result.
+ const real j = sign ? -OF : OF;
+ x = (j + x) - j;
+ exp = (vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff;
+
+ version (LittleEndian)
+ {
+ if (exp < 0)
+ result = 0;
+ else if (exp <= 31)
+ result = vi[1] >> (31 - exp);
+ else
+ result = (cast(long) vi[1] << (exp - 31)) | (vi[0] >> (63 - exp));
+ }
+ else
+ {
+ if (exp < 0)
+ result = 0;
+ else if (exp <= 31)
+ result = vi[1] >> (31 - exp);
+ else
+ result = (cast(long) vi[1] << (exp - 31)) | (vi[2] >> (63 - exp));
+ }
+ }
+ else
+ {
+ // It is left implementation defined when the number is too large
+ // to fit in a 64bit long.
+ return cast(long) x;
+ }
+
+ return sign ? -result : result;
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ const vu = cast(ushort*)(&x);
+
+ // Find the exponent and sign
+ const sign = (vu[F.EXPPOS_SHORT] >> 15) & 1;
+ if ((vu[F.EXPPOS_SHORT] & F.EXPMASK) - (F.EXPBIAS + 1) > 63)
+ {
+ // The result is left implementation defined when the number is
+ // too large to fit in a 64 bit long.
+ return cast(long) x;
+ }
+
+ // Force rounding of lower bits according to current rounding
+ // mode by adding ±2^-112 and subtracting it again.
+ enum OF = 5.19229685853482762853049632922009600E33L;
+ const j = sign ? -OF : OF;
+ x = (j + x) - j;
+
+ const implicitOne = 1UL << 48;
+ auto vl = cast(ulong*)(&x);
+ vl[MANTISSA_MSB] &= implicitOne - 1;
+ vl[MANTISSA_MSB] |= implicitOne;
+
+ long result;
+
+ const exp = (vu[F.EXPPOS_SHORT] & F.EXPMASK) - (F.EXPBIAS + 1);
+ if (exp < 0)
+ result = 0;
+ else if (exp <= 48)
+ result = vl[MANTISSA_MSB] >> (48 - exp);
+ else
+ result = (vl[MANTISSA_MSB] << (exp - 48)) | (vl[MANTISSA_LSB] >> (112 - exp));
+
+ return sign ? -result : result;
+ }
+ else
+ {
+ static assert(false, "real type not supported by lrint()");
+ }
+ }
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ assert(lrint(4.5) == 4);
+ assert(lrint(5.5) == 6);
+ assert(lrint(-4.5) == -4);
+ assert(lrint(-5.5) == -6);
+
+ assert(lrint(int.max - 0.5) == 2147483646L);
+ assert(lrint(int.max + 0.5) == 2147483648L);
+ assert(lrint(int.min - 0.5) == -2147483648L);
+ assert(lrint(int.min + 0.5) == -2147483648L);
+}
+
+static if (real.mant_dig >= long.sizeof * 8)
+{
+ @safe pure nothrow @nogc unittest
+ {
+ assert(lrint(long.max - 1.5L) == long.max - 1);
+ assert(lrint(long.max - 0.5L) == long.max - 1);
+ assert(lrint(long.min + 0.5L) == long.min);
+ assert(lrint(long.min + 1.5L) == long.min + 2);
+ }
+}
+
+/*******************************************
+ * Return the value of x rounded to the nearest integer.
+ * If the fractional part of x is exactly 0.5, the return value is
+ * rounded away from zero.
+ */
+real round(real x) @trusted nothrow @nogc
+{
+ version (CRuntime_Microsoft)
+ {
+ auto old = FloatingPointControl.getControlState();
+ FloatingPointControl.setControlState(
+ (old & ~FloatingPointControl.roundingMask) | FloatingPointControl.roundToZero
+ );
+ x = rint((x >= 0) ? x + 0.5 : x - 0.5);
+ FloatingPointControl.setControlState(old);
+ return x;
+ }
+ else
+ return core.stdc.math.roundl(x);
+}
+
+/**********************************************
+ * Return the value of x rounded to the nearest integer.
+ *
+ * If the fractional part of x is exactly 0.5, the return value is rounded
+ * away from zero.
+ *
+ * $(BLUE This function is Posix-Only.)
+ */
+long lround(real x) @trusted nothrow @nogc
+{
+ version (Posix)
+ return core.stdc.math.llroundl(x);
+ else
+ assert(0, "lround not implemented");
+}
+
+version (Posix)
+{
+ @safe nothrow @nogc unittest
+ {
+ assert(lround(0.49) == 0);
+ assert(lround(0.5) == 1);
+ assert(lround(1.5) == 2);
+ }
+}
+
+/****************************************************
+ * Returns the integer portion of x, dropping the fractional portion.
+ *
+ * This is also known as "chop" rounding.
+ */
+real trunc(real x) @trusted 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,0x0C ; // round to 0
+ 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,0x0C ; // round to 0
+ mov byte ptr cw+1,AL ;
+ fldcw cw ;
+ frndint ;
+ mov byte ptr cw+1,DL ;
+ fldcw cw ;
+ }
+ }
+ else
+ return core.stdc.math.truncl(x);
+}
+
+/****************************************************
+ * Calculate the remainder x REM y, following IEC 60559.
+ *
+ * REM is the value of x - y * n, where n is the integer nearest the exact
+ * value of x / y.
+ * If |n - x / y| == 0.5, n is even.
+ * If the result is zero, it has the same sign as x.
+ * Otherwise, the sign of the result is the sign of x / y.
+ * Precision mode has no effect on the remainder functions.
+ *
+ * remquo returns n in the parameter n.
+ *
+ * $(TABLE_SV
+ * $(TR $(TH x) $(TH y) $(TH remainder(x, y)) $(TH n) $(TH invalid?))
+ * $(TR $(TD $(PLUSMN)0.0) $(TD not 0.0) $(TD $(PLUSMN)0.0) $(TD 0.0) $(TD no))
+ * $(TR $(TD $(PLUSMNINF)) $(TD anything) $(TD $(NAN)) $(TD ?) $(TD yes))
+ * $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD $(NAN)) $(TD ?) $(TD yes))
+ * $(TR $(TD != $(PLUSMNINF)) $(TD $(PLUSMNINF)) $(TD x) $(TD ?) $(TD no))
+ * )
+ *
+ * $(BLUE `remquo` and `remainder` not supported on Windows.)
+ */
+real remainder(real x, real y) @trusted nothrow @nogc
+{
+ version (CRuntime_Microsoft)
+ {
+ int n;
+ return remquo(x, y, n);
+ }
+ else
+ return core.stdc.math.remainderl(x, y);
+}
+
+real remquo(real x, real y, out int n) @trusted nothrow @nogc /// ditto
+{
+ version (Posix)
+ return core.stdc.math.remquol(x, y, &n);
+ else
+ assert(0, "remquo not implemented");
+}
+
+/** IEEE exception status flags ('sticky bits')
+
+ These flags indicate that an exceptional floating-point condition has occurred.
+ They indicate that a NaN or an infinity has been generated, that a result
+ is inexact, or that a signalling NaN has been encountered. If floating-point
+ exceptions are enabled (unmasked), a hardware exception will be generated
+ instead of setting these flags.
+ */
+struct IeeeFlags
+{
+private:
+ // The x87 FPU status register is 16 bits.
+ // The Pentium SSE2 status register is 32 bits.
+ // The ARM and PowerPC FPSCR is a 32-bit register.
+ // The SPARC FSR is a 32bit register (64 bits for SPARC 7 & 8, but high bits are uninteresting).
+ uint flags;
+
+ version (CRuntime_Microsoft)
+ {
+ // Microsoft uses hardware-incompatible custom constants in fenv.h (core.stdc.fenv).
+ // Applies to both x87 status word (16 bits) and SSE2 status word(32 bits).
+ enum : int
+ {
+ INEXACT_MASK = 0x20,
+ UNDERFLOW_MASK = 0x10,
+ OVERFLOW_MASK = 0x08,
+ DIVBYZERO_MASK = 0x04,
+ INVALID_MASK = 0x01,
+
+ EXCEPTIONS_MASK = 0b11_1111
+ }
+ // Don't bother about subnormals, they are not supported on most CPUs.
+ // SUBNORMAL_MASK = 0x02;
+ }
+ else
+ {
+ enum : int
+ {
+ INEXACT_MASK = core.stdc.fenv.FE_INEXACT,
+ UNDERFLOW_MASK = core.stdc.fenv.FE_UNDERFLOW,
+ OVERFLOW_MASK = core.stdc.fenv.FE_OVERFLOW,
+ DIVBYZERO_MASK = core.stdc.fenv.FE_DIVBYZERO,
+ INVALID_MASK = core.stdc.fenv.FE_INVALID,
+ EXCEPTIONS_MASK = core.stdc.fenv.FE_ALL_EXCEPT,
+ }
+ }
+
+private:
+ static uint getIeeeFlags()
+ {
+ version (GNU)
+ {
+ version (X86_Any)
+ {
+ ushort sw;
+ asm pure nothrow @nogc
+ {
+ "fstsw %0" : "=a" (sw);
+ }
+ // OR the result with the SSE2 status register (MXCSR).
+ if (haveSSE)
+ {
+ uint mxcsr;
+ asm pure nothrow @nogc
+ {
+ "stmxcsr %0" : "=m" (mxcsr);
+ }
+ return (sw | mxcsr) & EXCEPTIONS_MASK;
+ }
+ else
+ return sw & EXCEPTIONS_MASK;
+ }
+ else version (ARM)
+ {
+ version (ARM_SoftFloat)
+ return 0;
+ else
+ {
+ uint result = void;
+ asm pure nothrow @nogc
+ {
+ "vmrs %0, FPSCR; and %0, %0, #0x1F;" : "=r" result;
+ }
+ return result;
+ }
+ }
+ else
+ assert(0, "Not yet supported");
+ }
+ else
+ version (InlineAsm_X86_Any)
+ {
+ ushort sw;
+ asm pure nothrow @nogc { fstsw sw; }
+
+ // OR the result with the SSE2 status register (MXCSR).
+ if (haveSSE)
+ {
+ uint mxcsr;
+ asm pure nothrow @nogc { stmxcsr mxcsr; }
+ return (sw | mxcsr) & EXCEPTIONS_MASK;
+ }
+ else return sw & EXCEPTIONS_MASK;
+ }
+ else version (SPARC)
+ {
+ /*
+ int retval;
+ asm pure nothrow @nogc { st %fsr, retval; }
+ return retval;
+ */
+ assert(0, "Not yet supported");
+ }
+ else version (ARM)
+ {
+ assert(false, "Not yet supported.");
+ }
+ else
+ assert(0, "Not yet supported");
+ }
+ static void resetIeeeFlags() @nogc
+ {
+ version (GNU)
+ {
+ version (X86_Any)
+ {
+ asm pure nothrow @nogc
+ {
+ "fnclex";
+ }
+
+ // Also clear exception flags in MXCSR, SSE's control register.
+ if (haveSSE)
+ {
+ uint mxcsr;
+ asm pure nothrow @nogc
+ {
+ "stmxcsr %0" : "=m" (mxcsr);
+ }
+ mxcsr &= ~EXCEPTIONS_MASK;
+ asm pure nothrow @nogc
+ {
+ "ldmxcsr %0" : : "m" (mxcsr);
+ }
+ }
+ }
+ else version (ARM)
+ {
+ version (ARM_SoftFloat)
+ return;
+ else
+ {
+ uint old = FloatingPointControl.getControlState();
+ old &= ~0b11111; // http://infocenter.arm.com/help/topic/com.arm.doc.ddi0408i/Chdfifdc.html
+ asm pure nothrow @nogc
+ {
+ "vmsr FPSCR, %0" : : "r" (old);
+ }
+ }
+ }
+ else
+ assert(0, "Not yet supported");
+ }
+ else
+ version (InlineAsm_X86_Any)
+ {
+ asm pure nothrow @nogc
+ {
+ fnclex;
+ }
+
+ // Also clear exception flags in MXCSR, SSE's control register.
+ if (haveSSE)
+ {
+ uint mxcsr;
+ asm nothrow @nogc { stmxcsr mxcsr; }
+ mxcsr &= ~EXCEPTIONS_MASK;
+ asm nothrow @nogc { ldmxcsr mxcsr; }
+ }
+ }
+ else
+ {
+ /* SPARC:
+ int tmpval;
+ asm pure nothrow @nogc { st %fsr, tmpval; }
+ tmpval &=0xFFFF_FC00;
+ asm pure nothrow @nogc { ld tmpval, %fsr; }
+ */
+ assert(0, "Not yet supported");
+ }
+ }
+public:
+ version (IeeeFlagsSupport)
+ {
+
+ /**
+ * The result cannot be represented exactly, so rounding occurred.
+ * Example: `x = sin(0.1);`
+ */
+ @property bool inexact() const { return (flags & INEXACT_MASK) != 0; }
+
+ /**
+ * A zero was generated by underflow
+ * Example: `x = real.min*real.epsilon/2;`
+ */
+ @property bool underflow() const { return (flags & UNDERFLOW_MASK) != 0; }
+
+ /**
+ * An infinity was generated by overflow
+ * Example: `x = real.max*2;`
+ */
+ @property bool overflow() const { return (flags & OVERFLOW_MASK) != 0; }
+
+ /**
+ * An infinity was generated by division by zero
+ * Example: `x = 3/0.0;`
+ */
+ @property bool divByZero() const { return (flags & DIVBYZERO_MASK) != 0; }
+
+ /**
+ * A machine NaN was generated.
+ * Example: `x = real.infinity * 0.0;`
+ */
+ @property bool invalid() const { return (flags & INVALID_MASK) != 0; }
+
+ }
+}
+
+///
+version (GNU)
+{
+ unittest
+ {
+ pragma(msg, "ieeeFlags test disabled, see LDC Issue #888");
+ }
+}
+else
+@system unittest
+{
+ static void func() {
+ int a = 10 * 10;
+ }
+
+ real a=3.5;
+ // Set all the flags to zero
+ resetIeeeFlags();
+ assert(!ieeeFlags.divByZero);
+ // Perform a division by zero.
+ a/=0.0L;
+ assert(a == real.infinity);
+ assert(ieeeFlags.divByZero);
+ // Create a NaN
+ a*=0.0L;
+ assert(ieeeFlags.invalid);
+ assert(isNaN(a));
+
+ // Check that calling func() has no effect on the
+ // status flags.
+ IeeeFlags f = ieeeFlags;
+ func();
+ assert(ieeeFlags == f);
+}
+
+version (GNU)
+{
+ unittest
+ {
+ pragma(msg, "ieeeFlags test disabled, see LDC Issue #888");
+ }
+}
+else
+@system unittest
+{
+ import std.meta : AliasSeq;
+
+ static struct Test
+ {
+ void delegate() action;
+ bool function() ieeeCheck;
+ }
+
+ foreach (T; AliasSeq!(float, double, real))
+ {
+ T x; /* Needs to be here to trick -O. It would optimize away the
+ calculations if x were local to the function literals. */
+ auto tests = [
+ Test(
+ () { x = 1; x += 0.1; },
+ () => ieeeFlags.inexact
+ ),
+ Test(
+ () { x = T.min_normal; x /= T.max; },
+ () => ieeeFlags.underflow
+ ),
+ Test(
+ () { x = T.max; x += T.max; },
+ () => ieeeFlags.overflow
+ ),
+ Test(
+ () { x = 1; x /= 0; },
+ () => ieeeFlags.divByZero
+ ),
+ Test(
+ () { x = 0; x /= 0; },
+ () => ieeeFlags.invalid
+ )
+ ];
+ foreach (test; tests)
+ {
+ resetIeeeFlags();
+ assert(!test.ieeeCheck());
+ test.action();
+ assert(test.ieeeCheck());
+ }
+ }
+}
+
+version (X86_Any)
+{
+ version = IeeeFlagsSupport;
+}
+version (X86_Any)
+{
+ version = IeeeFlagsSupport;
+}
+else version (PPC_Any)
+{
+ version = IeeeFlagsSupport;
+}
+else version (MIPS_Any)
+{
+ version = IeeeFlagsSupport;
+}
+else version (ARM_Any)
+{
+ version = IeeeFlagsSupport;
+}
+
+/// Set all of the floating-point status flags to false.
+void resetIeeeFlags() @nogc { IeeeFlags.resetIeeeFlags(); }
+
+/// Returns: snapshot of the current state of the floating-point status flags
+@property IeeeFlags ieeeFlags()
+{
+ return IeeeFlags(IeeeFlags.getIeeeFlags());
+}
+
+/** Control the Floating point hardware
+
+ Change the IEEE754 floating-point rounding mode and the floating-point
+ hardware exceptions.
+
+ By default, the rounding mode is roundToNearest and all hardware exceptions
+ are disabled. For most applications, debugging is easier if the $(I division
+ by zero), $(I overflow), and $(I invalid operation) exceptions are enabled.
+ These three are combined into a $(I severeExceptions) value for convenience.
+ Note in particular that if $(I invalidException) is enabled, a hardware trap
+ will be generated whenever an uninitialized floating-point variable is used.
+
+ All changes are temporary. The previous state is restored at the
+ end of the scope.
+
+
+Example:
+----
+{
+ FloatingPointControl fpctrl;
+
+ // Enable hardware exceptions for division by zero, overflow to infinity,
+ // invalid operations, and uninitialized floating-point variables.
+ fpctrl.enableExceptions(FloatingPointControl.severeExceptions);
+
+ // This will generate a hardware exception, if x is a
+ // default-initialized floating point variable:
+ real x; // Add `= 0` or even `= real.nan` to not throw the exception.
+ real y = x * 3.0;
+
+ // The exception is only thrown for default-uninitialized NaN-s.
+ // NaN-s with other payload are valid:
+ real z = y * real.nan; // ok
+
+ // Changing the rounding mode:
+ fpctrl.rounding = FloatingPointControl.roundUp;
+ assert(rint(1.1) == 2);
+
+ // The set hardware exceptions will be disabled when leaving this scope.
+ // The original rounding mode will also be restored.
+}
+
+// Ensure previous values are returned:
+assert(!FloatingPointControl.enabledExceptions);
+assert(FloatingPointControl.rounding == FloatingPointControl.roundToNearest);
+assert(rint(1.1) == 1);
+----
+
+ */
+struct FloatingPointControl
+{
+ alias RoundingMode = uint; ///
+
+ version (StdDdoc)
+ {
+ enum : RoundingMode
+ {
+ /** IEEE rounding modes.
+ * The default mode is roundToNearest.
+ *
+ * roundingMask = A mask of all rounding modes.
+ */
+ roundToNearest,
+ roundDown, /// ditto
+ roundUp, /// ditto
+ roundToZero, /// ditto
+ roundingMask, /// ditto
+ }
+ }
+ else version (CRuntime_Microsoft)
+ {
+ // Microsoft uses hardware-incompatible custom constants in fenv.h (core.stdc.fenv).
+ enum : RoundingMode
+ {
+ roundToNearest = 0x0000,
+ roundDown = 0x0400,
+ roundUp = 0x0800,
+ roundToZero = 0x0C00,
+ roundingMask = roundToNearest | roundDown
+ | roundUp | roundToZero,
+ }
+ }
+ else
+ {
+ enum : RoundingMode
+ {
+ roundToNearest = core.stdc.fenv.FE_TONEAREST,
+ roundDown = core.stdc.fenv.FE_DOWNWARD,
+ roundUp = core.stdc.fenv.FE_UPWARD,
+ roundToZero = core.stdc.fenv.FE_TOWARDZERO,
+ roundingMask = roundToNearest | roundDown
+ | roundUp | roundToZero,
+ }
+ }
+
+ //// Change the floating-point hardware rounding mode
+ @property void rounding(RoundingMode newMode) @nogc
+ {
+ initialize();
+ setControlState(cast(ushort)((getControlState() & (-1 - roundingMask)) | (newMode & roundingMask)));
+ }
+
+ /// Returns: the currently active rounding mode
+ @property static RoundingMode rounding() @nogc
+ {
+ return cast(RoundingMode)(getControlState() & roundingMask);
+ }
+
+ alias ExceptionMask = uint; ///
+
+ version (StdDdoc)
+ {
+ enum : ExceptionMask
+ {
+ /** IEEE hardware exceptions.
+ * By default, all exceptions are masked (disabled).
+ *
+ * severeExceptions = The overflow, division by zero, and invalid
+ * exceptions.
+ */
+ subnormalException,
+ inexactException, /// ditto
+ underflowException, /// ditto
+ overflowException, /// ditto
+ divByZeroException, /// ditto
+ invalidException, /// ditto
+ severeExceptions, /// ditto
+ allExceptions, /// ditto
+ }
+ }
+ else version (ARM_Any)
+ {
+ enum : ExceptionMask
+ {
+ subnormalException = 0x8000,
+ inexactException = 0x1000,
+ underflowException = 0x0800,
+ overflowException = 0x0400,
+ divByZeroException = 0x0200,
+ invalidException = 0x0100,
+ severeExceptions = overflowException | divByZeroException
+ | invalidException,
+ allExceptions = severeExceptions | underflowException
+ | inexactException | subnormalException,
+ }
+ }
+ else version (MIPS_Any)
+ {
+ enum : ExceptionMask
+ {
+ inexactException = 0x0080,
+ underflowException = 0x0100,
+ overflowException = 0x0200,
+ divByZeroException = 0x0400,
+ invalidException = 0x0800,
+ severeExceptions = overflowException | divByZeroException
+ | invalidException,
+ allExceptions = severeExceptions | underflowException
+ | inexactException,
+ }
+ }
+ else version (PPC_Any)
+ {
+ enum : ExceptionMask
+ {
+ inexactException = 0x08,
+ divByZeroException = 0x10,
+ underflowException = 0x20,
+ overflowException = 0x40,
+ invalidException = 0x80,
+ severeExceptions = overflowException | divByZeroException
+ | invalidException,
+ allExceptions = severeExceptions | underflowException
+ | inexactException,
+ }
+ }
+ else version (SPARC64)
+ {
+ enum : ExceptionMask
+ {
+ inexactException = 0x0800000,
+ divByZeroException = 0x1000000,
+ overflowException = 0x4000000,
+ underflowException = 0x2000000,
+ invalidException = 0x8000000,
+ severeExceptions = overflowException | divByZeroException
+ | invalidException,
+ allExceptions = severeExceptions | underflowException
+ | inexactException,
+ }
+ }
+ else version (SystemZ)
+ {
+ enum : ExceptionMask
+ {
+ inexactException = 0x08000000,
+ divByZeroException = 0x40000000,
+ overflowException = 0x20000000,
+ underflowException = 0x10000000,
+ invalidException = 0x80000000,
+ severeExceptions = overflowException | divByZeroException
+ | invalidException,
+ allExceptions = severeExceptions | underflowException
+ | inexactException,
+ }
+ }
+ else version (X86_Any)
+ {
+ enum : ExceptionMask
+ {
+ inexactException = 0x20,
+ underflowException = 0x10,
+ overflowException = 0x08,
+ divByZeroException = 0x04,
+ subnormalException = 0x02,
+ invalidException = 0x01,
+ severeExceptions = overflowException | divByZeroException
+ | invalidException,
+ allExceptions = severeExceptions | underflowException
+ | inexactException | subnormalException,
+ }
+ }
+ else
+ static assert(false, "Not implemented for this architecture");
+
+public:
+ /// Returns: true if the current FPU supports exception trapping
+ @property static bool hasExceptionTraps() @safe nothrow @nogc
+ {
+ version (X86_Any)
+ return true;
+ else version (PPC_Any)
+ return true;
+ else version (MIPS_Any)
+ return true;
+ else version (ARM_Any)
+ {
+ auto oldState = getControlState();
+ // If exceptions are not supported, we set the bit but read it back as zero
+ // https://sourceware.org/ml/libc-ports/2012-06/msg00091.html
+ setControlState(oldState | divByZeroException);
+ immutable result = (getControlState() & allExceptions) != 0;
+ setControlState(oldState);
+ return result;
+ }
+ else
+ assert(0, "Not yet supported");
+ }
+
+ /// Enable (unmask) specific hardware exceptions. Multiple exceptions may be ORed together.
+ void enableExceptions(ExceptionMask exceptions) @nogc
+ {
+ assert(hasExceptionTraps);
+ initialize();
+ version (X86_Any)
+ setControlState(getControlState() & ~(exceptions & allExceptions));
+ else
+ setControlState(getControlState() | (exceptions & allExceptions));
+ }
+
+ /// Disable (mask) specific hardware exceptions. Multiple exceptions may be ORed together.
+ void disableExceptions(ExceptionMask exceptions) @nogc
+ {
+ assert(hasExceptionTraps);
+ initialize();
+ version (X86_Any)
+ setControlState(getControlState() | (exceptions & allExceptions));
+ else
+ setControlState(getControlState() & ~(exceptions & allExceptions));
+ }
+
+ /// Returns: the exceptions which are currently enabled (unmasked)
+ @property static ExceptionMask enabledExceptions() @nogc
+ {
+ assert(hasExceptionTraps);
+ version (X86_Any)
+ return (getControlState() & allExceptions) ^ allExceptions;
+ else
+ return (getControlState() & allExceptions);
+ }
+
+ /// Clear all pending exceptions, then restore the original exception state and rounding mode.
+ ~this() @nogc
+ {
+ clearExceptions();
+ if (initialized)
+ setControlState(savedState);
+ }
+
+private:
+ ControlState savedState;
+
+ bool initialized = false;
+
+ version (ARM_Any)
+ {
+ alias ControlState = uint;
+ }
+ else version (PPC_Any)
+ {
+ alias ControlState = uint;
+ }
+ else version (MIPS_Any)
+ {
+ alias ControlState = uint;
+ }
+ else version (SPARC64)
+ {
+ alias ControlState = ulong;
+ }
+ else version (SystemZ)
+ {
+ alias ControlState = uint;
+ }
+ else version (X86_Any)
+ {
+ alias ControlState = ushort;
+ }
+ else
+ static assert(false, "Not implemented for this architecture");
+
+ void initialize() @nogc
+ {
+ // BUG: This works around the absence of this() constructors.
+ if (initialized) return;
+ clearExceptions();
+ savedState = getControlState();
+ initialized = true;
+ }
+
+ // Clear all pending exceptions
+ static void clearExceptions() @nogc
+ {
+ resetIeeeFlags();
+ }
+
+ // Read from the control register
+ static ControlState getControlState() @trusted nothrow @nogc
+ {
+ version (GNU)
+ {
+ version (X86_Any)
+ {
+ ControlState cont;
+ asm pure nothrow @nogc
+ {
+ "fstcw %0" : "=m" cont;
+ }
+ return cont;
+ }
+ else version (AArch64)
+ {
+ asm pure nothrow @nogc
+ {
+ "mrs %0, FPCR;" : "=r" cont;
+ }
+ return cont;
+ }
+ else version (ARM)
+ {
+ ControlState cont;
+ version (ARM_SoftFloat)
+ cont = 0;
+ else
+ {
+ asm pure nothrow @nogc
+ {
+ "vmrs %0, FPSCR" : "=r" cont;
+ }
+ }
+ return cont;
+ }
+ else
+ assert(0, "Not yet supported");
+ }
+ else
+ version (D_InlineAsm_X86)
+ {
+ short cont;
+ asm nothrow @nogc
+ {
+ xor EAX, EAX;
+ fstcw cont;
+ }
+ return cont;
+ }
+ else
+ version (D_InlineAsm_X86_64)
+ {
+ short cont;
+ asm nothrow @nogc
+ {
+ xor RAX, RAX;
+ fstcw cont;
+ }
+ return cont;
+ }
+ else
+ assert(0, "Not yet supported");
+ }
+
+ // Set the control register
+ static void setControlState(ControlState newState) @trusted nothrow @nogc
+ {
+ version (GNU)
+ {
+ version (X86_Any)
+ {
+ asm pure nothrow @nogc
+ {
+ "fclex; fldcw %0" : : "m" newState;
+ }
+
+ // Also update MXCSR, SSE's control register.
+ if (haveSSE)
+ {
+ uint mxcsr;
+ asm pure nothrow @nogc
+ {
+ "stmxcsr %0" : "=m" mxcsr;
+ }
+
+ /* In the FPU control register, rounding mode is in bits 10 and
+ 11. In MXCSR it's in bits 13 and 14. */
+ mxcsr &= ~(roundingMask << 3); // delete old rounding mode
+ mxcsr |= (newState & roundingMask) << 3; // write new rounding mode
+
+ /* In the FPU control register, masks are bits 0 through 5.
+ In MXCSR they're 7 through 12. */
+ mxcsr &= ~(allExceptions << 7); // delete old masks
+ mxcsr |= (newState & allExceptions) << 7; // write new exception masks
+
+ asm pure nothrow @nogc
+ {
+ "ldmxcsr %0" : : "m" mxcsr;
+ }
+ }
+ }
+ else version (AArch64)
+ {
+ asm pure nothrow @nogc
+ {
+ "msr FPCR, %0;" : : "r" (newState);
+ }
+ }
+ else version (ARM)
+ {
+ version (ARM_SoftFloat)
+ return;
+ else
+ {
+ asm pure nothrow @nogc
+ {
+ "vmsr FPSCR, %0" : : "r" (newState);
+ }
+ }
+ }
+ else
+ assert(0, "Not yet supported");
+ }
+ else
+ version (InlineAsm_X86_Any)
+ {
+ asm nothrow @nogc
+ {
+ fclex;
+ fldcw newState;
+ }
+
+ // Also update MXCSR, SSE's control register.
+ if (haveSSE)
+ {
+ uint mxcsr;
+ asm nothrow @nogc { stmxcsr mxcsr; }
+
+ /* In the FPU control register, rounding mode is in bits 10 and
+ 11. In MXCSR it's in bits 13 and 14. */
+ mxcsr &= ~(roundingMask << 3); // delete old rounding mode
+ mxcsr |= (newState & roundingMask) << 3; // write new rounding mode
+
+ /* In the FPU control register, masks are bits 0 through 5.
+ In MXCSR they're 7 through 12. */
+ mxcsr &= ~(allExceptions << 7); // delete old masks
+ mxcsr |= (newState & allExceptions) << 7; // write new exception masks
+
+ asm nothrow @nogc { ldmxcsr mxcsr; }
+ }
+ }
+ else
+ assert(0, "Not yet supported");
+ }
+}
+
+@system unittest
+{
+ // GCC floating point emulation doesn't allow changing
+ // rounding modes, getting error bits etc
+ version (GNU) version (D_SoftFloat)
+ return;
+
+ void ensureDefaults()
+ {
+ assert(FloatingPointControl.rounding
+ == FloatingPointControl.roundToNearest);
+ if (FloatingPointControl.hasExceptionTraps)
+ assert(FloatingPointControl.enabledExceptions == 0);
+ }
+
+ {
+ FloatingPointControl ctrl;
+ }
+ ensureDefaults();
+
+ version (D_HardFloat)
+ {
+ {
+ FloatingPointControl ctrl;
+ ctrl.rounding = FloatingPointControl.roundDown;
+ assert(FloatingPointControl.rounding == FloatingPointControl.roundDown);
+ }
+ ensureDefaults();
+ }
+
+ if (FloatingPointControl.hasExceptionTraps)
+ {
+ FloatingPointControl ctrl;
+ ctrl.enableExceptions(FloatingPointControl.divByZeroException
+ | FloatingPointControl.overflowException);
+ assert(ctrl.enabledExceptions ==
+ (FloatingPointControl.divByZeroException
+ | FloatingPointControl.overflowException));
+
+ ctrl.rounding = FloatingPointControl.roundUp;
+ assert(FloatingPointControl.rounding == FloatingPointControl.roundUp);
+ }
+ ensureDefaults();
+}
+
+@system unittest // rounding
+{
+ import std.meta : AliasSeq;
+
+ foreach (T; AliasSeq!(float, double, real))
+ {
+ FloatingPointControl fpctrl;
+
+ fpctrl.rounding = FloatingPointControl.roundUp;
+ T u = 1;
+ u += 0.1;
+
+ fpctrl.rounding = FloatingPointControl.roundDown;
+ T d = 1;
+ d += 0.1;
+
+ fpctrl.rounding = FloatingPointControl.roundToZero;
+ T z = 1;
+ z += 0.1;
+
+ assert(u > d);
+ assert(z == d);
+
+ fpctrl.rounding = FloatingPointControl.roundUp;
+ u = -1;
+ u -= 0.1;
+
+ fpctrl.rounding = FloatingPointControl.roundDown;
+ d = -1;
+ d -= 0.1;
+
+ fpctrl.rounding = FloatingPointControl.roundToZero;
+ z = -1;
+ z -= 0.1;
+
+ assert(u > d);
+ assert(z == u);
+ }
+}
+
+
+/*********************************
+ * Determines if $(D_PARAM x) is NaN.
+ * Params:
+ * x = a floating point number.
+ * Returns:
+ * $(D true) if $(D_PARAM x) is Nan.
+ */
+bool isNaN(X)(X x) @nogc @trusted pure nothrow
+if (isFloatingPoint!(X))
+{
+ alias F = floatTraits!(X);
+ static if (F.realFormat == RealFormat.ieeeSingle)
+ {
+ const uint p = *cast(uint *)&x;
+ return ((p & 0x7F80_0000) == 0x7F80_0000)
+ && p & 0x007F_FFFF; // not infinity
+ }
+ else static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ const ulong p = *cast(ulong *)&x;
+ return ((p & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000)
+ && p & 0x000F_FFFF_FFFF_FFFF; // not infinity
+ }
+ else static if (F.realFormat == RealFormat.ieeeExtended)
+ {
+ const ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
+ const ulong ps = *cast(ulong *)&x;
+ return e == F.EXPMASK &&
+ ps & 0x7FFF_FFFF_FFFF_FFFF; // not infinity
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ const ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
+ const ulong psLsb = (cast(ulong *)&x)[MANTISSA_LSB];
+ const ulong psMsb = (cast(ulong *)&x)[MANTISSA_MSB];
+ return e == F.EXPMASK &&
+ (psLsb | (psMsb& 0x0000_FFFF_FFFF_FFFF)) != 0;
+ }
+ else
+ {
+ return x != x;
+ }
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ assert( isNaN(float.init));
+ assert( isNaN(-double.init));
+ assert( isNaN(real.nan));
+ assert( isNaN(-real.nan));
+ assert(!isNaN(cast(float) 53.6));
+ assert(!isNaN(cast(real)-53.6));
+}
+
+@safe pure nothrow @nogc unittest
+{
+ import std.meta : AliasSeq;
+
+ foreach (T; AliasSeq!(float, double, real))
+ {
+ // CTFE-able tests
+ assert(isNaN(T.init));
+ assert(isNaN(-T.init));
+ assert(isNaN(T.nan));
+ assert(isNaN(-T.nan));
+ assert(!isNaN(T.infinity));
+ assert(!isNaN(-T.infinity));
+ assert(!isNaN(cast(T) 53.6));
+ assert(!isNaN(cast(T)-53.6));
+
+ // Runtime tests
+ shared T f;
+ f = T.init;
+ assert(isNaN(f));
+ assert(isNaN(-f));
+ f = T.nan;
+ assert(isNaN(f));
+ assert(isNaN(-f));
+ f = T.infinity;
+ assert(!isNaN(f));
+ assert(!isNaN(-f));
+ f = cast(T) 53.6;
+ assert(!isNaN(f));
+ assert(!isNaN(-f));
+ }
+}
+
+/*********************************
+ * Determines if $(D_PARAM x) is finite.
+ * Params:
+ * x = a floating point number.
+ * Returns:
+ * $(D true) if $(D_PARAM x) is finite.
+ */
+bool isFinite(X)(X x) @trusted pure nothrow @nogc
+{
+ alias F = floatTraits!(X);
+ ushort* pe = cast(ushort *)&x;
+ return (pe[F.EXPPOS_SHORT] & F.EXPMASK) != F.EXPMASK;
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ assert( isFinite(1.23f));
+ assert( isFinite(float.max));
+ assert( isFinite(float.min_normal));
+ assert(!isFinite(float.nan));
+ assert(!isFinite(float.infinity));
+}
+
+@safe pure nothrow @nogc unittest
+{
+ assert(isFinite(1.23));
+ assert(isFinite(double.max));
+ assert(isFinite(double.min_normal));
+ assert(!isFinite(double.nan));
+ assert(!isFinite(double.infinity));
+
+ assert(isFinite(1.23L));
+ assert(isFinite(real.max));
+ assert(isFinite(real.min_normal));
+ assert(!isFinite(real.nan));
+ assert(!isFinite(real.infinity));
+}
+
+
+/*********************************
+ * Determines if $(D_PARAM x) is normalized.
+ *
+ * A normalized number must not be zero, subnormal, infinite nor $(NAN).
+ *
+ * Params:
+ * x = a floating point number.
+ * Returns:
+ * $(D true) if $(D_PARAM x) is normalized.
+ */
+
+/* Need one for each format because subnormal floats might
+ * be converted to normal reals.
+ */
+bool isNormal(X)(X x) @trusted pure nothrow @nogc
+{
+ alias F = floatTraits!(X);
+ static if (F.realFormat == RealFormat.ibmExtended)
+ {
+ // doubledouble is normal if the least significant part is normal.
+ return isNormal((cast(double*)&x)[MANTISSA_LSB]);
+ }
+ else
+ {
+ ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
+ return (e != F.EXPMASK && e != 0);
+ }
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ float f = 3;
+ double d = 500;
+ real e = 10e+48;
+
+ assert(isNormal(f));
+ assert(isNormal(d));
+ assert(isNormal(e));
+ f = d = e = 0;
+ assert(!isNormal(f));
+ assert(!isNormal(d));
+ assert(!isNormal(e));
+ assert(!isNormal(real.infinity));
+ assert(isNormal(-real.max));
+ assert(!isNormal(real.min_normal/4));
+
+}
+
+/*********************************
+ * Determines if $(D_PARAM x) is subnormal.
+ *
+ * Subnormals (also known as "denormal number"), have a 0 exponent
+ * and a 0 most significant mantissa bit.
+ *
+ * Params:
+ * x = a floating point number.
+ * Returns:
+ * $(D true) if $(D_PARAM x) is a denormal number.
+ */
+bool isSubnormal(X)(X x) @trusted pure nothrow @nogc
+{
+ /*
+ Need one for each format because subnormal floats might
+ be converted to normal reals.
+ */
+ alias F = floatTraits!(X);
+ static if (F.realFormat == RealFormat.ieeeSingle)
+ {
+ uint *p = cast(uint *)&x;
+ return (*p & F.EXPMASK_INT) == 0 && *p & F.MANTISSAMASK_INT;
+ }
+ else static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ uint *p = cast(uint *)&x;
+ return (p[MANTISSA_MSB] & F.EXPMASK_INT) == 0
+ && (p[MANTISSA_LSB] || p[MANTISSA_MSB] & F.MANTISSAMASK_INT);
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
+ long* ps = cast(long *)&x;
+ return (e == 0 &&
+ ((ps[MANTISSA_LSB]|(ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF)) != 0));
+ }
+ else static if (F.realFormat == RealFormat.ieeeExtended)
+ {
+ ushort* pe = cast(ushort *)&x;
+ long* ps = cast(long *)&x;
+
+ return (pe[F.EXPPOS_SHORT] & F.EXPMASK) == 0 && *ps > 0;
+ }
+ else static if (F.realFormat == RealFormat.ibmExtended)
+ {
+ return isSubnormal((cast(double*)&x)[MANTISSA_MSB]);
+ }
+ else
+ {
+ static assert(false, "Not implemented for this architecture");
+ }
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ import std.meta : AliasSeq;
+
+ foreach (T; AliasSeq!(float, double, real))
+ {
+ T f;
+ for (f = 1.0; !isSubnormal(f); f /= 2)
+ assert(f != 0);
+ }
+}
+
+/*********************************
+ * Determines if $(D_PARAM x) is $(PLUSMN)$(INFIN).
+ * Params:
+ * x = a floating point number.
+ * Returns:
+ * $(D true) if $(D_PARAM x) is $(PLUSMN)$(INFIN).
+ */
+bool isInfinity(X)(X x) @nogc @trusted pure nothrow
+if (isFloatingPoint!(X))
+{
+ alias F = floatTraits!(X);
+ static if (F.realFormat == RealFormat.ieeeSingle)
+ {
+ return ((*cast(uint *)&x) & 0x7FFF_FFFF) == 0x7F80_0000;
+ }
+ else static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF)
+ == 0x7FF0_0000_0000_0000;
+ }
+ else static if (F.realFormat == RealFormat.ieeeExtended)
+ {
+ const ushort e = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]);
+ const ulong ps = *cast(ulong *)&x;
+
+ // On Motorola 68K, infinity can have hidden bit = 1 or 0. On x86, it is always 1.
+ return e == F.EXPMASK && (ps & 0x7FFF_FFFF_FFFF_FFFF) == 0;
+ }
+ else static if (F.realFormat == RealFormat.ibmExtended)
+ {
+ return (((cast(ulong *)&x)[MANTISSA_MSB]) & 0x7FFF_FFFF_FFFF_FFFF)
+ == 0x7FF8_0000_0000_0000;
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ const long psLsb = (cast(long *)&x)[MANTISSA_LSB];
+ const long psMsb = (cast(long *)&x)[MANTISSA_MSB];
+ return (psLsb == 0)
+ && (psMsb & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_0000_0000_0000;
+ }
+ else
+ {
+ return (x < -X.max) || (X.max < x);
+ }
+}
+
+///
+@nogc @safe pure nothrow unittest
+{
+ assert(!isInfinity(float.init));
+ assert(!isInfinity(-float.init));
+ assert(!isInfinity(float.nan));
+ assert(!isInfinity(-float.nan));
+ assert(isInfinity(float.infinity));
+ assert(isInfinity(-float.infinity));
+ assert(isInfinity(-1.0f / 0.0f));
+}
+
+@safe pure nothrow @nogc unittest
+{
+ // CTFE-able tests
+ assert(!isInfinity(double.init));
+ assert(!isInfinity(-double.init));
+ assert(!isInfinity(double.nan));
+ assert(!isInfinity(-double.nan));
+ assert(isInfinity(double.infinity));
+ assert(isInfinity(-double.infinity));
+ assert(isInfinity(-1.0 / 0.0));
+
+ assert(!isInfinity(real.init));
+ assert(!isInfinity(-real.init));
+ assert(!isInfinity(real.nan));
+ assert(!isInfinity(-real.nan));
+ assert(isInfinity(real.infinity));
+ assert(isInfinity(-real.infinity));
+ assert(isInfinity(-1.0L / 0.0L));
+
+ // Runtime tests
+ shared float f;
+ f = float.init;
+ assert(!isInfinity(f));
+ assert(!isInfinity(-f));
+ f = float.nan;
+ assert(!isInfinity(f));
+ assert(!isInfinity(-f));
+ f = float.infinity;
+ assert(isInfinity(f));
+ assert(isInfinity(-f));
+ f = (-1.0f / 0.0f);
+ assert(isInfinity(f));
+
+ shared double d;
+ d = double.init;
+ assert(!isInfinity(d));
+ assert(!isInfinity(-d));
+ d = double.nan;
+ assert(!isInfinity(d));
+ assert(!isInfinity(-d));
+ d = double.infinity;
+ assert(isInfinity(d));
+ assert(isInfinity(-d));
+ d = (-1.0 / 0.0);
+ assert(isInfinity(d));
+
+ shared real e;
+ e = real.init;
+ assert(!isInfinity(e));
+ assert(!isInfinity(-e));
+ e = real.nan;
+ assert(!isInfinity(e));
+ assert(!isInfinity(-e));
+ e = real.infinity;
+ assert(isInfinity(e));
+ assert(isInfinity(-e));
+ e = (-1.0L / 0.0L);
+ assert(isInfinity(e));
+}
+
+/*********************************
+ * Is the binary representation of x identical to y?
+ *
+ * Same as ==, except that positive and negative zero are not identical,
+ * and two $(NAN)s are identical if they have the same 'payload'.
+ */
+bool isIdentical(real x, real y) @trusted pure nothrow @nogc
+{
+ // We're doing a bitwise comparison so the endianness is irrelevant.
+ long* pxs = cast(long *)&x;
+ long* pys = cast(long *)&y;
+ alias F = floatTraits!(real);
+ static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ return pxs[0] == pys[0];
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple
+ || F.realFormat == RealFormat.ibmExtended)
+ {
+ return pxs[0] == pys[0] && pxs[1] == pys[1];
+ }
+ else
+ {
+ ushort* pxe = cast(ushort *)&x;
+ ushort* pye = cast(ushort *)&y;
+ return pxe[4] == pye[4] && pxs[0] == pys[0];
+ }
+}
+
+/*********************************
+ * Return 1 if sign bit of e is set, 0 if not.
+ */
+int signbit(X)(X x) @nogc @trusted pure nothrow
+{
+ alias F = floatTraits!(X);
+ return ((cast(ubyte *)&x)[F.SIGNPOS_BYTE] & 0x80) != 0;
+}
+
+///
+@nogc @safe pure nothrow unittest
+{
+ assert(!signbit(float.nan));
+ assert(signbit(-float.nan));
+ assert(!signbit(168.1234f));
+ assert(signbit(-168.1234f));
+ assert(!signbit(0.0f));
+ assert(signbit(-0.0f));
+ assert(signbit(-float.max));
+ assert(!signbit(float.max));
+
+ assert(!signbit(double.nan));
+ assert(signbit(-double.nan));
+ assert(!signbit(168.1234));
+ assert(signbit(-168.1234));
+ assert(!signbit(0.0));
+ assert(signbit(-0.0));
+ assert(signbit(-double.max));
+ assert(!signbit(double.max));
+
+ assert(!signbit(real.nan));
+ assert(signbit(-real.nan));
+ assert(!signbit(168.1234L));
+ assert(signbit(-168.1234L));
+ assert(!signbit(0.0L));
+ assert(signbit(-0.0L));
+ assert(signbit(-real.max));
+ assert(!signbit(real.max));
+}
+
+
+/*********************************
+ * Return a value composed of to with from's sign bit.
+ */
+R copysign(R, X)(R to, X from) @trusted pure nothrow @nogc
+if (isFloatingPoint!(R) && isFloatingPoint!(X))
+{
+ ubyte* pto = cast(ubyte *)&to;
+ const ubyte* pfrom = cast(ubyte *)&from;
+
+ alias T = floatTraits!(R);
+ alias F = floatTraits!(X);
+ pto[T.SIGNPOS_BYTE] &= 0x7F;
+ pto[T.SIGNPOS_BYTE] |= pfrom[F.SIGNPOS_BYTE] & 0x80;
+ return to;
+}
+
+// ditto
+R copysign(R, X)(X to, R from) @trusted pure nothrow @nogc
+if (isIntegral!(X) && isFloatingPoint!(R))
+{
+ return copysign(cast(R) to, from);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ import std.meta : AliasSeq;
+
+ foreach (X; AliasSeq!(float, double, real, int, long))
+ {
+ foreach (Y; AliasSeq!(float, double, real))
+ (){ // avoid slow optimizations for large functions @@@BUG@@@ 2396
+ X x = 21;
+ Y y = 23.8;
+ Y e = void;
+
+ e = copysign(x, y);
+ assert(e == 21.0);
+
+ e = copysign(-x, y);
+ assert(e == 21.0);
+
+ e = copysign(x, -y);
+ assert(e == -21.0);
+
+ e = copysign(-x, -y);
+ assert(e == -21.0);
+
+ static if (isFloatingPoint!X)
+ {
+ e = copysign(X.nan, y);
+ assert(isNaN(e) && !signbit(e));
+
+ e = copysign(X.nan, -y);
+ assert(isNaN(e) && signbit(e));
+ }
+ }();
+ }
+}
+
+/*********************************
+Returns $(D -1) if $(D x < 0), $(D x) if $(D x == 0), $(D 1) if
+$(D x > 0), and $(NAN) if x==$(NAN).
+ */
+F sgn(F)(F x) @safe pure nothrow @nogc
+{
+ // @@@TODO@@@: make this faster
+ return x > 0 ? 1 : x < 0 ? -1 : x;
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ assert(sgn(168.1234) == 1);
+ assert(sgn(-168.1234) == -1);
+ assert(sgn(0.0) == 0);
+ assert(sgn(-0.0) == 0);
+}
+
+// Functions for NaN payloads
+/*
+ * A 'payload' can be stored in the significand of a $(NAN). One bit is required
+ * to distinguish between a quiet and a signalling $(NAN). This leaves 22 bits
+ * of payload for a float; 51 bits for a double; 62 bits for an 80-bit real;
+ * and 111 bits for a 128-bit quad.
+*/
+/**
+ * Create a quiet $(NAN), storing an integer inside the payload.
+ *
+ * For floats, the largest possible payload is 0x3F_FFFF.
+ * For doubles, it is 0x3_FFFF_FFFF_FFFF.
+ * For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF.
+ */
+real NaN(ulong payload) @trusted pure nothrow @nogc
+{
+ alias F = floatTraits!(real);
+ static if (F.realFormat == RealFormat.ieeeExtended)
+ {
+ // real80 (in x86 real format, the implied bit is actually
+ // not implied but a real bit which is stored in the real)
+ ulong v = 3; // implied bit = 1, quiet bit = 1
+ }
+ else
+ {
+ ulong v = 1; // no implied bit. quiet bit = 1
+ }
+
+ ulong a = payload;
+
+ // 22 Float bits
+ ulong w = a & 0x3F_FFFF;
+ a -= w;
+
+ v <<=22;
+ v |= w;
+ a >>=22;
+
+ // 29 Double bits
+ v <<=29;
+ w = a & 0xFFF_FFFF;
+ v |= w;
+ a -= w;
+ a >>=29;
+
+ static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ v |= 0x7FF0_0000_0000_0000;
+ real x;
+ * cast(ulong *)(&x) = v;
+ return x;
+ }
+ else
+ {
+ v <<=11;
+ a &= 0x7FF;
+ v |= a;
+ real x = real.nan;
+
+ // Extended real bits
+ static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ v <<= 1; // there's no implicit bit
+
+ version (LittleEndian)
+ {
+ *cast(ulong*)(6+cast(ubyte*)(&x)) = v;
+ }
+ else
+ {
+ *cast(ulong*)(2+cast(ubyte*)(&x)) = v;
+ }
+ }
+ else
+ {
+ *cast(ulong *)(&x) = v;
+ }
+ return x;
+ }
+}
+
+@system pure nothrow @nogc unittest // not @safe because taking address of local.
+{
+ static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
+ {
+ auto x = NaN(1);
+ auto xl = *cast(ulong*)&x;
+ assert(xl & 0x8_0000_0000_0000UL); //non-signaling bit, bit 52
+ assert((xl & 0x7FF0_0000_0000_0000UL) == 0x7FF0_0000_0000_0000UL); //all exp bits set
+ }
+}
+
+/**
+ * Extract an integral payload from a $(NAN).
+ *
+ * Returns:
+ * the integer payload as a ulong.
+ *
+ * For floats, the largest possible payload is 0x3F_FFFF.
+ * For doubles, it is 0x3_FFFF_FFFF_FFFF.
+ * For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF.
+ */
+ulong getNaNPayload(real x) @trusted pure nothrow @nogc
+{
+ // assert(isNaN(x));
+ alias F = floatTraits!(real);
+ static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ ulong m = *cast(ulong *)(&x);
+ // Make it look like an 80-bit significand.
+ // Skip exponent, and quiet bit
+ m &= 0x0007_FFFF_FFFF_FFFF;
+ m <<= 11;
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ version (LittleEndian)
+ {
+ ulong m = *cast(ulong*)(6+cast(ubyte*)(&x));
+ }
+ else
+ {
+ ulong m = *cast(ulong*)(2+cast(ubyte*)(&x));
+ }
+
+ m >>= 1; // there's no implicit bit
+ }
+ else
+ {
+ ulong m = *cast(ulong *)(&x);
+ }
+
+ // ignore implicit bit and quiet bit
+
+ const ulong f = m & 0x3FFF_FF00_0000_0000L;
+
+ ulong w = f >>> 40;
+ w |= (m & 0x00FF_FFFF_F800L) << (22 - 11);
+ w |= (m & 0x7FF) << 51;
+ return w;
+}
+
+debug(UnitTest)
+{
+ @safe pure nothrow @nogc unittest
+ {
+ real nan4 = NaN(0x789_ABCD_EF12_3456);
+ static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended
+ || floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
+ {
+ assert(getNaNPayload(nan4) == 0x789_ABCD_EF12_3456);
+ }
+ else
+ {
+ assert(getNaNPayload(nan4) == 0x1_ABCD_EF12_3456);
+ }
+ double nan5 = nan4;
+ assert(getNaNPayload(nan5) == 0x1_ABCD_EF12_3456);
+ float nan6 = nan4;
+ assert(getNaNPayload(nan6) == 0x12_3456);
+ nan4 = NaN(0xFABCD);
+ assert(getNaNPayload(nan4) == 0xFABCD);
+ nan6 = nan4;
+ assert(getNaNPayload(nan6) == 0xFABCD);
+ nan5 = NaN(0x100_0000_0000_3456);
+ assert(getNaNPayload(nan5) == 0x0000_0000_3456);
+ }
+}
+
+/**
+ * Calculate the next largest floating point value after x.
+ *
+ * Return the least number greater than x that is representable as a real;
+ * thus, it gives the next point on the IEEE number line.
+ *
+ * $(TABLE_SV
+ * $(SVH x, nextUp(x) )
+ * $(SV -$(INFIN), -real.max )
+ * $(SV $(PLUSMN)0.0, real.min_normal*real.epsilon )
+ * $(SV real.max, $(INFIN) )
+ * $(SV $(INFIN), $(INFIN) )
+ * $(SV $(NAN), $(NAN) )
+ * )
+ */
+real nextUp(real x) @trusted pure nothrow @nogc
+{
+ alias F = floatTraits!(real);
+ static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ return nextUp(cast(double) x);
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
+ if (e == F.EXPMASK)
+ {
+ // NaN or Infinity
+ if (x == -real.infinity) return -real.max;
+ return x; // +Inf and NaN are unchanged.
+ }
+
+ auto ps = cast(ulong *)&x;
+ if (ps[MANTISSA_MSB] & 0x8000_0000_0000_0000)
+ {
+ // Negative number
+ if (ps[MANTISSA_LSB] == 0 && ps[MANTISSA_MSB] == 0x8000_0000_0000_0000)
+ {
+ // it was negative zero, change to smallest subnormal
+ ps[MANTISSA_LSB] = 1;
+ ps[MANTISSA_MSB] = 0;
+ return x;
+ }
+ if (ps[MANTISSA_LSB] == 0) --ps[MANTISSA_MSB];
+ --ps[MANTISSA_LSB];
+ }
+ else
+ {
+ // Positive number
+ ++ps[MANTISSA_LSB];
+ if (ps[MANTISSA_LSB] == 0) ++ps[MANTISSA_MSB];
+ }
+ return x;
+ }
+ else static if (F.realFormat == RealFormat.ieeeExtended)
+ {
+ // For 80-bit reals, the "implied bit" is a nuisance...
+ ushort *pe = cast(ushort *)&x;
+ ulong *ps = cast(ulong *)&x;
+
+ if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK)
+ {
+ // First, deal with NANs and infinity
+ if (x == -real.infinity) return -real.max;
+ return x; // +Inf and NaN are unchanged.
+ }
+ if (pe[F.EXPPOS_SHORT] & 0x8000)
+ {
+ // Negative number -- need to decrease the significand
+ --*ps;
+ // Need to mask with 0x7FFF... so subnormals are treated correctly.
+ if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF)
+ {
+ if (pe[F.EXPPOS_SHORT] == 0x8000) // it was negative zero
+ {
+ *ps = 1;
+ pe[F.EXPPOS_SHORT] = 0; // smallest subnormal.
+ return x;
+ }
+
+ --pe[F.EXPPOS_SHORT];
+
+ if (pe[F.EXPPOS_SHORT] == 0x8000)
+ return x; // it's become a subnormal, implied bit stays low.
+
+ *ps = 0xFFFF_FFFF_FFFF_FFFF; // set the implied bit
+ return x;
+ }
+ return x;
+ }
+ else
+ {
+ // Positive number -- need to increase the significand.
+ // Works automatically for positive zero.
+ ++*ps;
+ if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0)
+ {
+ // change in exponent
+ ++pe[F.EXPPOS_SHORT];
+ *ps = 0x8000_0000_0000_0000; // set the high bit
+ }
+ }
+ return x;
+ }
+ else // static if (F.realFormat == RealFormat.ibmExtended)
+ {
+ assert(0, "nextUp not implemented");
+ }
+}
+
+/** ditto */
+double nextUp(double x) @trusted pure nothrow @nogc
+{
+ ulong *ps = cast(ulong *)&x;
+
+ if ((*ps & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000)
+ {
+ // First, deal with NANs and infinity
+ if (x == -x.infinity) return -x.max;
+ return x; // +INF and NAN are unchanged.
+ }
+ if (*ps & 0x8000_0000_0000_0000) // Negative number
+ {
+ if (*ps == 0x8000_0000_0000_0000) // it was negative zero
+ {
+ *ps = 0x0000_0000_0000_0001; // change to smallest subnormal
+ return x;
+ }
+ --*ps;
+ }
+ else
+ { // Positive number
+ ++*ps;
+ }
+ return x;
+}
+
+/** ditto */
+float nextUp(float x) @trusted pure nothrow @nogc
+{
+ uint *ps = cast(uint *)&x;
+
+ if ((*ps & 0x7F80_0000) == 0x7F80_0000)
+ {
+ // First, deal with NANs and infinity
+ if (x == -x.infinity) return -x.max;
+
+ return x; // +INF and NAN are unchanged.
+ }
+ if (*ps & 0x8000_0000) // Negative number
+ {
+ if (*ps == 0x8000_0000) // it was negative zero
+ {
+ *ps = 0x0000_0001; // change to smallest subnormal
+ return x;
+ }
+
+ --*ps;
+ }
+ else
+ {
+ // Positive number
+ ++*ps;
+ }
+ return x;
+}
+
+/**
+ * Calculate the next smallest floating point value before x.
+ *
+ * Return the greatest number less than x that is representable as a real;
+ * thus, it gives the previous point on the IEEE number line.
+ *
+ * $(TABLE_SV
+ * $(SVH x, nextDown(x) )
+ * $(SV $(INFIN), real.max )
+ * $(SV $(PLUSMN)0.0, -real.min_normal*real.epsilon )
+ * $(SV -real.max, -$(INFIN) )
+ * $(SV -$(INFIN), -$(INFIN) )
+ * $(SV $(NAN), $(NAN) )
+ * )
+ */
+real nextDown(real x) @safe pure nothrow @nogc
+{
+ return -nextUp(-x);
+}
+
+/** ditto */
+double nextDown(double x) @safe pure nothrow @nogc
+{
+ return -nextUp(-x);
+}
+
+/** ditto */
+float nextDown(float x) @safe pure nothrow @nogc
+{
+ return -nextUp(-x);
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ assert( nextDown(1.0 + real.epsilon) == 1.0);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
+ {
+
+ // Tests for 80-bit reals
+ assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC)));
+ // negative numbers
+ assert( nextUp(-real.infinity) == -real.max );
+ assert( nextUp(-1.0L-real.epsilon) == -1.0 );
+ assert( nextUp(-2.0L) == -2.0 + real.epsilon);
+ // subnormals and zero
+ assert( nextUp(-real.min_normal) == -real.min_normal*(1-real.epsilon) );
+ assert( nextUp(-real.min_normal*(1-real.epsilon)) == -real.min_normal*(1-2*real.epsilon) );
+ assert( isIdentical(-0.0L, nextUp(-real.min_normal*real.epsilon)) );
+ assert( nextUp(-0.0L) == real.min_normal*real.epsilon );
+ assert( nextUp(0.0L) == real.min_normal*real.epsilon );
+ assert( nextUp(real.min_normal*(1-real.epsilon)) == real.min_normal );
+ assert( nextUp(real.min_normal) == real.min_normal*(1+real.epsilon) );
+ // positive numbers
+ assert( nextUp(1.0L) == 1.0 + real.epsilon );
+ assert( nextUp(2.0L-real.epsilon) == 2.0 );
+ assert( nextUp(real.max) == real.infinity );
+ assert( nextUp(real.infinity)==real.infinity );
+ }
+
+ double n = NaN(0xABC);
+ assert(isIdentical(nextUp(n), n));
+ // negative numbers
+ assert( nextUp(-double.infinity) == -double.max );
+ assert( nextUp(-1-double.epsilon) == -1.0 );
+ assert( nextUp(-2.0) == -2.0 + double.epsilon);
+ // subnormals and zero
+
+ assert( nextUp(-double.min_normal) == -double.min_normal*(1-double.epsilon) );
+ assert( nextUp(-double.min_normal*(1-double.epsilon)) == -double.min_normal*(1-2*double.epsilon) );
+ assert( isIdentical(-0.0, nextUp(-double.min_normal*double.epsilon)) );
+ assert( nextUp(0.0) == double.min_normal*double.epsilon );
+ assert( nextUp(-0.0) == double.min_normal*double.epsilon );
+ assert( nextUp(double.min_normal*(1-double.epsilon)) == double.min_normal );
+ assert( nextUp(double.min_normal) == double.min_normal*(1+double.epsilon) );
+ // positive numbers
+ assert( nextUp(1.0) == 1.0 + double.epsilon );
+ assert( nextUp(2.0-double.epsilon) == 2.0 );
+ assert( nextUp(double.max) == double.infinity );
+
+ float fn = NaN(0xABC);
+ assert(isIdentical(nextUp(fn), fn));
+ float f = -float.min_normal*(1-float.epsilon);
+ float f1 = -float.min_normal;
+ assert( nextUp(f1) == f);
+ f = 1.0f+float.epsilon;
+ f1 = 1.0f;
+ assert( nextUp(f1) == f );
+ f1 = -0.0f;
+ assert( nextUp(f1) == float.min_normal*float.epsilon);
+ assert( nextUp(float.infinity)==float.infinity );
+
+ assert(nextDown(1.0L+real.epsilon)==1.0);
+ assert(nextDown(1.0+double.epsilon)==1.0);
+ f = 1.0f+float.epsilon;
+ assert(nextDown(f)==1.0);
+ assert(nextafter(1.0+real.epsilon, -real.infinity)==1.0);
+}
+
+
+
+/******************************************
+ * Calculates the next representable value after x in the direction of y.
+ *
+ * If y > x, the result will be the next largest floating-point value;
+ * if y < x, the result will be the next smallest value.
+ * If x == y, the result is y.
+ *
+ * Remarks:
+ * This function is not generally very useful; it's almost always better to use
+ * the faster functions nextUp() or nextDown() instead.
+ *
+ * The FE_INEXACT and FE_OVERFLOW exceptions will be raised if x is finite and
+ * the function result is infinite. The FE_INEXACT and FE_UNDERFLOW
+ * exceptions will be raised if the function value is subnormal, and x is
+ * not equal to y.
+ */
+T nextafter(T)(const T x, const T y) @safe pure nothrow @nogc
+{
+ if (x == y) return y;
+ return ((y>x) ? nextUp(x) : nextDown(x));
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ float a = 1;
+ assert(is(typeof(nextafter(a, a)) == float));
+ assert(nextafter(a, a.infinity) > a);
+
+ double b = 2;
+ assert(is(typeof(nextafter(b, b)) == double));
+ assert(nextafter(b, b.infinity) > b);
+
+ real c = 3;
+ assert(is(typeof(nextafter(c, c)) == real));
+ assert(nextafter(c, c.infinity) > c);
+}
+
+//real nexttoward(real x, real y) { return core.stdc.math.nexttowardl(x, y); }
+
+/*******************************************
+ * Returns the positive difference between x and y.
+ * Returns:
+ * $(TABLE_SV
+ * $(TR $(TH x, y) $(TH fdim(x, y)))
+ * $(TR $(TD x $(GT) y) $(TD x - y))
+ * $(TR $(TD x $(LT)= y) $(TD +0.0))
+ * )
+ */
+real fdim(real x, real y) @safe pure nothrow @nogc { return (x > y) ? x - y : +0.0; }
+
+/****************************************
+ * Returns the larger of x and y.
+ */
+real fmax(real x, real y) @safe pure nothrow @nogc { return x > y ? x : y; }
+
+/****************************************
+ * Returns the smaller of x and y.
+ */
+real fmin(real x, real y) @safe pure nothrow @nogc { return x < y ? x : y; }
+
+/**************************************
+ * Returns (x * y) + z, rounding only once according to the
+ * current rounding mode.
+ *
+ * BUGS: Not currently implemented - rounds twice.
+ */
+real fma(real x, real y, real z) @safe pure nothrow @nogc { return (x * y) + z; }
+
+/*******************************************************************
+ * Compute the value of x $(SUPERSCRIPT n), where n is an integer
+ */
+Unqual!F pow(F, G)(F x, G n) @nogc @trusted pure nothrow
+if (isFloatingPoint!(F) && isIntegral!(G))
+{
+ import std.traits : Unsigned;
+ real p = 1.0, v = void;
+ Unsigned!(Unqual!G) m = n;
+ if (n < 0)
+ {
+ switch (n)
+ {
+ case -1:
+ return 1 / x;
+ case -2:
+ return 1 / (x * x);
+ default:
+ }
+
+ m = cast(typeof(m))(0 - n);
+ v = p / x;
+ }
+ else
+ {
+ switch (n)
+ {
+ case 0:
+ return 1.0;
+ case 1:
+ return x;
+ case 2:
+ return x * x;
+ default:
+ }
+
+ v = x;
+ }
+
+ while (1)
+ {
+ if (m & 1)
+ p *= v;
+ m >>= 1;
+ if (!m)
+ break;
+ v *= v;
+ }
+ return p;
+}
+
+@safe pure nothrow @nogc unittest
+{
+ // Make sure it instantiates and works properly on immutable values and
+ // with various integer and float types.
+ immutable real x = 46;
+ immutable float xf = x;
+ immutable double xd = x;
+ immutable uint one = 1;
+ immutable ushort two = 2;
+ immutable ubyte three = 3;
+ immutable ulong eight = 8;
+
+ immutable int neg1 = -1;
+ immutable short neg2 = -2;
+ immutable byte neg3 = -3;
+ immutable long neg8 = -8;
+
+
+ assert(pow(x,0) == 1.0);
+ assert(pow(xd,one) == x);
+ assert(pow(xf,two) == x * x);
+ assert(pow(x,three) == x * x * x);
+ assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x));
+
+ assert(pow(x, neg1) == 1 / x);
+
+ version (X86_64)
+ {
+ pragma(msg, "test disabled on x86_64, see bug 5628");
+ }
+ else version (ARM)
+ {
+ pragma(msg, "test disabled on ARM, see bug 5628");
+ }
+ else
+ {
+ assert(pow(xd, neg2) == 1 / (x * x));
+ assert(pow(xf, neg8) == 1 / ((x * x) * (x * x) * (x * x) * (x * x)));
+ }
+
+ assert(feqrel(pow(x, neg3), 1 / (x * x * x)) >= real.mant_dig - 1);
+}
+
+@system unittest
+{
+ assert(equalsDigit(pow(2.0L, 10.0L), 1024, 19));
+}
+
+/** Compute the value of an integer x, raised to the power of a positive
+ * integer n.
+ *
+ * If both x and n are 0, the result is 1.
+ * If n is negative, an integer divide error will occur at runtime,
+ * regardless of the value of x.
+ */
+typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @nogc @trusted pure nothrow
+if (isIntegral!(F) && isIntegral!(G))
+{
+ if (n<0) return x/0; // Only support positive powers
+ typeof(return) p, v = void;
+ Unqual!G m = n;
+
+ switch (m)
+ {
+ case 0:
+ p = 1;
+ break;
+
+ case 1:
+ p = x;
+ break;
+
+ case 2:
+ p = x * x;
+ break;
+
+ default:
+ v = x;
+ p = 1;
+ while (1)
+ {
+ if (m & 1)
+ p *= v;
+ m >>= 1;
+ if (!m)
+ break;
+ v *= v;
+ }
+ break;
+ }
+ return p;
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ immutable int one = 1;
+ immutable byte two = 2;
+ immutable ubyte three = 3;
+ immutable short four = 4;
+ immutable long ten = 10;
+
+ assert(pow(two, three) == 8);
+ assert(pow(two, ten) == 1024);
+ assert(pow(one, ten) == 1);
+ assert(pow(ten, four) == 10_000);
+ assert(pow(four, 10) == 1_048_576);
+ assert(pow(three, four) == 81);
+
+}
+
+/**Computes integer to floating point powers.*/
+real pow(I, F)(I x, F y) @nogc @trusted pure nothrow
+if (isIntegral!I && isFloatingPoint!F)
+{
+ return pow(cast(real) x, cast(Unqual!F) y);
+}
+
+/*********************************************
+ * Calculates x$(SUPERSCRIPT y).
+ *
+ * $(TABLE_SV
+ * $(TR $(TH x) $(TH y) $(TH pow(x, y))
+ * $(TH div 0) $(TH invalid?))
+ * $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD 1.0)
+ * $(TD no) $(TD no) )
+ * $(TR $(TD |x| $(GT) 1) $(TD +$(INFIN)) $(TD +$(INFIN))
+ * $(TD no) $(TD no) )
+ * $(TR $(TD |x| $(LT) 1) $(TD +$(INFIN)) $(TD +0.0)
+ * $(TD no) $(TD no) )
+ * $(TR $(TD |x| $(GT) 1) $(TD -$(INFIN)) $(TD +0.0)
+ * $(TD no) $(TD no) )
+ * $(TR $(TD |x| $(LT) 1) $(TD -$(INFIN)) $(TD +$(INFIN))
+ * $(TD no) $(TD no) )
+ * $(TR $(TD +$(INFIN)) $(TD $(GT) 0.0) $(TD +$(INFIN))
+ * $(TD no) $(TD no) )
+ * $(TR $(TD +$(INFIN)) $(TD $(LT) 0.0) $(TD +0.0)
+ * $(TD no) $(TD no) )
+ * $(TR $(TD -$(INFIN)) $(TD odd integer $(GT) 0.0) $(TD -$(INFIN))
+ * $(TD no) $(TD no) )
+ * $(TR $(TD -$(INFIN)) $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN))
+ * $(TD no) $(TD no))
+ * $(TR $(TD -$(INFIN)) $(TD odd integer $(LT) 0.0) $(TD -0.0)
+ * $(TD no) $(TD no) )
+ * $(TR $(TD -$(INFIN)) $(TD $(LT) 0.0, not odd integer) $(TD +0.0)
+ * $(TD no) $(TD no) )
+ * $(TR $(TD $(PLUSMN)1.0) $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))
+ * $(TD no) $(TD yes) )
+ * $(TR $(TD $(LT) 0.0) $(TD finite, nonintegral) $(TD $(NAN))
+ * $(TD no) $(TD yes))
+ * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(LT) 0.0) $(TD $(PLUSMNINF))
+ * $(TD yes) $(TD no) )
+ * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN))
+ * $(TD yes) $(TD no))
+ * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(GT) 0.0) $(TD $(PLUSMN)0.0)
+ * $(TD no) $(TD no) )
+ * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT) 0.0, not odd integer) $(TD +0.0)
+ * $(TD no) $(TD no) )
+ * )
+ */
+Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow
+if (isFloatingPoint!(F) && isFloatingPoint!(G))
+{
+ alias Float = typeof(return);
+
+ static real impl(real x, real y) @nogc pure nothrow
+ {
+ // Special cases.
+ if (isNaN(y))
+ return y;
+ if (isNaN(x) && y != 0.0)
+ return x;
+
+ // Even if x is NaN.
+ if (y == 0.0)
+ return 1.0;
+ if (y == 1.0)
+ return x;
+
+ if (isInfinity(y))
+ {
+ if (fabs(x) > 1)
+ {
+ if (signbit(y))
+ return +0.0;
+ else
+ return F.infinity;
+ }
+ else if (fabs(x) == 1)
+ {
+ return y * 0; // generate NaN.
+ }
+ else // < 1
+ {
+ if (signbit(y))
+ return F.infinity;
+ else
+ return +0.0;
+ }
+ }
+ if (isInfinity(x))
+ {
+ if (signbit(x))
+ {
+ long i = cast(long) y;
+ if (y > 0.0)
+ {
+ if (i == y && i & 1)
+ return -F.infinity;
+ else
+ return F.infinity;
+ }
+ else if (y < 0.0)
+ {
+ if (i == y && i & 1)
+ return -0.0;
+ else
+ return +0.0;
+ }
+ }
+ else
+ {
+ if (y > 0.0)
+ return F.infinity;
+ else if (y < 0.0)
+ return +0.0;
+ }
+ }
+
+ if (x == 0.0)
+ {
+ if (signbit(x))
+ {
+ long i = cast(long) y;
+ if (y > 0.0)
+ {
+ if (i == y && i & 1)
+ return -0.0;
+ else
+ return +0.0;
+ }
+ else if (y < 0.0)
+ {
+ if (i == y && i & 1)
+ return -F.infinity;
+ else
+ return F.infinity;
+ }
+ }
+ else
+ {
+ if (y > 0.0)
+ return +0.0;
+ else if (y < 0.0)
+ return F.infinity;
+ }
+ }
+ if (x == 1.0)
+ return 1.0;
+
+ if (y >= F.max)
+ {
+ if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
+ return 0.0;
+ if (x > 1.0 || x < -1.0)
+ return F.infinity;
+ }
+ if (y <= -F.max)
+ {
+ if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
+ return F.infinity;
+ if (x > 1.0 || x < -1.0)
+ return 0.0;
+ }
+
+ if (x >= F.max)
+ {
+ if (y > 0.0)
+ return F.infinity;
+ else
+ return 0.0;
+ }
+ if (x <= -F.max)
+ {
+ long i = cast(long) y;
+ if (y > 0.0)
+ {
+ if (i == y && i & 1)
+ return -F.infinity;
+ else
+ return F.infinity;
+ }
+ else if (y < 0.0)
+ {
+ if (i == y && i & 1)
+ return -0.0;
+ else
+ return +0.0;
+ }
+ }
+
+ // Integer power of x.
+ long iy = cast(long) y;
+ if (iy == y && fabs(y) < 32_768.0)
+ return pow(x, iy);
+
+ real sign = 1.0;
+ if (x < 0)
+ {
+ // Result is real only if y is an integer
+ // Check for a non-zero fractional part
+ enum maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
+ static if (maxOdd > ulong.max)
+ {
+ // Generic method, for any FP type
+ if (floor(y) != y)
+ return sqrt(x); // Complex result -- create a NaN
+
+ const hy = ldexp(y, -1);
+ if (floor(hy) != hy)
+ sign = -1.0;
+ }
+ else
+ {
+ // Much faster, if ulong has enough precision
+ const absY = fabs(y);
+ if (absY <= maxOdd)
+ {
+ const uy = cast(ulong) absY;
+ if (uy != absY)
+ return sqrt(x); // Complex result -- create a NaN
+
+ if (uy & 1)
+ sign = -1.0;
+ }
+ }
+ x = -x;
+ }
+ version (INLINE_YL2X)
+ {
+ // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
+ // TODO: This is not accurate in practice. A fast and accurate
+ // (though complicated) method is described in:
+ // "An efficient rounding boundary test for pow(x, y)
+ // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
+ return sign * exp2( core.math.yl2x(x, y) );
+ }
+ else
+ {
+ // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
+ // TODO: This is not accurate in practice. A fast and accurate
+ // (though complicated) method is described in:
+ // "An efficient rounding boundary test for pow(x, y)
+ // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
+ Float w = exp2(y * log2(x));
+ return sign * w;
+ }
+ }
+ return impl(x, y);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ // Test all the special values. These unittests can be run on Windows
+ // by temporarily changing the version (linux) to version (all).
+ immutable float zero = 0;
+ immutable real one = 1;
+ immutable double two = 2;
+ immutable float three = 3;
+ immutable float fnan = float.nan;
+ immutable double dnan = double.nan;
+ immutable real rnan = real.nan;
+ immutable dinf = double.infinity;
+ immutable rninf = -real.infinity;
+
+ assert(pow(fnan, zero) == 1);
+ assert(pow(dnan, zero) == 1);
+ assert(pow(rnan, zero) == 1);
+
+ assert(pow(two, dinf) == double.infinity);
+ assert(isIdentical(pow(0.2f, dinf), +0.0));
+ assert(pow(0.99999999L, rninf) == real.infinity);
+ assert(isIdentical(pow(1.000000001, rninf), +0.0));
+ assert(pow(dinf, 0.001) == dinf);
+ assert(isIdentical(pow(dinf, -0.001), +0.0));
+ assert(pow(rninf, 3.0L) == rninf);
+ assert(pow(rninf, 2.0L) == real.infinity);
+ assert(isIdentical(pow(rninf, -3.0), -0.0));
+ assert(isIdentical(pow(rninf, -2.0), +0.0));
+
+ // @@@BUG@@@ somewhere
+ version (OSX) {} else assert(isNaN(pow(one, dinf)));
+ version (OSX) {} else assert(isNaN(pow(-one, dinf)));
+ assert(isNaN(pow(-0.2, PI)));
+ // boundary cases. Note that epsilon == 2^^-n for some n,
+ // so 1/epsilon == 2^^n is always even.
+ assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L);
+ assert(pow(-1.0L, 1/real.epsilon) == 1.0L);
+ assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L)));
+ assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L)));
+
+ assert(pow(0.0, -3.0) == double.infinity);
+ assert(pow(-0.0, -3.0) == -double.infinity);
+ assert(pow(0.0, -PI) == double.infinity);
+ assert(pow(-0.0, -PI) == double.infinity);
+ assert(isIdentical(pow(0.0, 5.0), 0.0));
+ assert(isIdentical(pow(-0.0, 5.0), -0.0));
+ assert(isIdentical(pow(0.0, 6.0), 0.0));
+ assert(isIdentical(pow(-0.0, 6.0), 0.0));
+
+ // Issue #14786 fixed
+ immutable real maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
+ assert(pow(-1.0L, maxOdd) == -1.0L);
+ assert(pow(-1.0L, -maxOdd) == -1.0L);
+ assert(pow(-1.0L, maxOdd + 1.0L) == 1.0L);
+ assert(pow(-1.0L, -maxOdd + 1.0L) == 1.0L);
+ assert(pow(-1.0L, maxOdd - 1.0L) == 1.0L);
+ assert(pow(-1.0L, -maxOdd - 1.0L) == 1.0L);
+
+ // Now, actual numbers.
+ assert(approxEqual(pow(two, three), 8.0));
+ assert(approxEqual(pow(two, -2.5), 0.1767767));
+
+ // Test integer to float power.
+ immutable uint twoI = 2;
+ assert(approxEqual(pow(twoI, three), 8.0));
+}
+
+/**************************************
+ * To what precision is x equal to y?
+ *
+ * Returns: the number of mantissa bits which are equal in x and y.
+ * eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision.
+ *
+ * $(TABLE_SV
+ * $(TR $(TH x) $(TH y) $(TH feqrel(x, y)))
+ * $(TR $(TD x) $(TD x) $(TD real.mant_dig))
+ * $(TR $(TD x) $(TD $(GT)= 2*x) $(TD 0))
+ * $(TR $(TD x) $(TD $(LT)= x/2) $(TD 0))
+ * $(TR $(TD $(NAN)) $(TD any) $(TD 0))
+ * $(TR $(TD any) $(TD $(NAN)) $(TD 0))
+ * )
+ */
+int feqrel(X)(const X x, const X y) @trusted pure nothrow @nogc
+if (isFloatingPoint!(X))
+{
+ /* Public Domain. Author: Don Clugston, 18 Aug 2005.
+ */
+ alias F = floatTraits!(X);
+ static if (F.realFormat == RealFormat.ibmExtended)
+ {
+ if (cast(double*)(&x)[MANTISSA_MSB] == cast(double*)(&y)[MANTISSA_MSB])
+ {
+ return double.mant_dig
+ + feqrel(cast(double*)(&x)[MANTISSA_LSB],
+ cast(double*)(&y)[MANTISSA_LSB]);
+ }
+ else
+ {
+ return feqrel(cast(double*)(&x)[MANTISSA_MSB],
+ cast(double*)(&y)[MANTISSA_MSB]);
+ }
+ }
+ else
+ {
+ static assert(F.realFormat == RealFormat.ieeeSingle
+ || F.realFormat == RealFormat.ieeeDouble
+ || F.realFormat == RealFormat.ieeeExtended
+ || F.realFormat == RealFormat.ieeeQuadruple);
+
+ if (x == y)
+ return X.mant_dig; // ensure diff != 0, cope with INF.
+
+ Unqual!X diff = fabs(x - y);
+
+ ushort *pa = cast(ushort *)(&x);
+ ushort *pb = cast(ushort *)(&y);
+ ushort *pd = cast(ushort *)(&diff);
+
+
+ // The difference in abs(exponent) between x or y and abs(x-y)
+ // is equal to the number of significand bits of x which are
+ // equal to y. If negative, x and y have different exponents.
+ // If positive, x and y are equal to 'bitsdiff' bits.
+ // AND with 0x7FFF to form the absolute value.
+ // To avoid out-by-1 errors, we subtract 1 so it rounds down
+ // if the exponents were different. This means 'bitsdiff' is
+ // always 1 lower than we want, except that if bitsdiff == 0,
+ // they could have 0 or 1 bits in common.
+
+ int bitsdiff = ((( (pa[F.EXPPOS_SHORT] & F.EXPMASK)
+ + (pb[F.EXPPOS_SHORT] & F.EXPMASK)
+ - (1 << F.EXPSHIFT)) >> 1)
+ - (pd[F.EXPPOS_SHORT] & F.EXPMASK)) >> F.EXPSHIFT;
+ if ( (pd[F.EXPPOS_SHORT] & F.EXPMASK) == 0)
+ { // Difference is subnormal
+ // For subnormals, we need to add the number of zeros that
+ // lie at the start of diff's significand.
+ // We do this by multiplying by 2^^real.mant_dig
+ diff *= F.RECIP_EPSILON;
+ return bitsdiff + X.mant_dig - ((pd[F.EXPPOS_SHORT] & F.EXPMASK) >> F.EXPSHIFT);
+ }
+
+ if (bitsdiff > 0)
+ return bitsdiff + 1; // add the 1 we subtracted before
+
+ // Avoid out-by-1 errors when factor is almost 2.
+ if (bitsdiff == 0
+ && ((pa[F.EXPPOS_SHORT] ^ pb[F.EXPPOS_SHORT]) & F.EXPMASK) == 0)
+ {
+ return 1;
+ } else return 0;
+ }
+}
+
+@safe pure nothrow @nogc unittest
+{
+ void testFeqrel(F)()
+ {
+ // Exact equality
+ assert(feqrel(F.max, F.max) == F.mant_dig);
+ assert(feqrel!(F)(0.0, 0.0) == F.mant_dig);
+ assert(feqrel(F.infinity, F.infinity) == F.mant_dig);
+
+ // a few bits away from exact equality
+ F w=1;
+ for (int i = 1; i < F.mant_dig - 1; ++i)
+ {
+ assert(feqrel!(F)(1.0 + w * F.epsilon, 1.0) == F.mant_dig-i);
+ assert(feqrel!(F)(1.0 - w * F.epsilon, 1.0) == F.mant_dig-i);
+ assert(feqrel!(F)(1.0, 1 + (w-1) * F.epsilon) == F.mant_dig - i + 1);
+ w*=2;
+ }
+
+ assert(feqrel!(F)(1.5+F.epsilon, 1.5) == F.mant_dig-1);
+ assert(feqrel!(F)(1.5-F.epsilon, 1.5) == F.mant_dig-1);
+ assert(feqrel!(F)(1.5-F.epsilon, 1.5+F.epsilon) == F.mant_dig-2);
+
+
+ // Numbers that are close
+ assert(feqrel!(F)(0x1.Bp+84, 0x1.B8p+84) == 5);
+ assert(feqrel!(F)(0x1.8p+10, 0x1.Cp+10) == 2);
+ assert(feqrel!(F)(1.5 * (1 - F.epsilon), 1.0L) == 2);
+ assert(feqrel!(F)(1.5, 1.0) == 1);
+ assert(feqrel!(F)(2 * (1 - F.epsilon), 1.0L) == 1);
+
+ // Factors of 2
+ assert(feqrel(F.max, F.infinity) == 0);
+ assert(feqrel!(F)(2 * (1 - F.epsilon), 1.0L) == 1);
+ assert(feqrel!(F)(1.0, 2.0) == 0);
+ assert(feqrel!(F)(4.0, 1.0) == 0);
+
+ // Extreme inequality
+ assert(feqrel(F.nan, F.nan) == 0);
+ assert(feqrel!(F)(0.0L, -F.nan) == 0);
+ assert(feqrel(F.nan, F.infinity) == 0);
+ assert(feqrel(F.infinity, -F.infinity) == 0);
+ assert(feqrel(F.max, -F.max) == 0);
+
+ assert(feqrel(F.min_normal / 8, F.min_normal / 17) == 3);
+
+ const F Const = 2;
+ immutable F Immutable = 2;
+ auto Compiles = feqrel(Const, Immutable);
+ }
+
+ assert(feqrel(7.1824L, 7.1824L) == real.mant_dig);
+
+ testFeqrel!(real)();
+ testFeqrel!(double)();
+ testFeqrel!(float)();
+}
+
+package: // Not public yet
+/* Return the value that lies halfway between x and y on the IEEE number line.
+ *
+ * Formally, the result is the arithmetic mean of the binary significands of x
+ * and y, multiplied by the geometric mean of the binary exponents of x and y.
+ * x and y must have the same sign, and must not be NaN.
+ * Note: this function is useful for ensuring O(log n) behaviour in algorithms
+ * involving a 'binary chop'.
+ *
+ * Special cases:
+ * If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value
+ * is the arithmetic mean (x + y) / 2.
+ * If x and y are even powers of 2, the return value is the geometric mean,
+ * ieeeMean(x, y) = sqrt(x * y).
+ *
+ */
+T ieeeMean(T)(const T x, const T y) @trusted pure nothrow @nogc
+in
+{
+ // both x and y must have the same sign, and must not be NaN.
+ assert(signbit(x) == signbit(y));
+ assert(x == x && y == y);
+}
+body
+{
+ // Runtime behaviour for contract violation:
+ // If signs are opposite, or one is a NaN, return 0.
+ if (!((x >= 0 && y >= 0) || (x <= 0 && y <= 0))) return 0.0;
+
+ // The implementation is simple: cast x and y to integers,
+ // average them (avoiding overflow), and cast the result back to a floating-point number.
+
+ alias F = floatTraits!(T);
+ T u;
+ static if (F.realFormat == RealFormat.ieeeExtended)
+ {
+ // There's slight additional complexity because they are actually
+ // 79-bit reals...
+ ushort *ue = cast(ushort *)&u;
+ ulong *ul = cast(ulong *)&u;
+ ushort *xe = cast(ushort *)&x;
+ ulong *xl = cast(ulong *)&x;
+ ushort *ye = cast(ushort *)&y;
+ ulong *yl = cast(ulong *)&y;
+
+ // Ignore the useless implicit bit. (Bonus: this prevents overflows)
+ ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL);
+
+ // @@@ BUG? @@@
+ // Cast shouldn't be here
+ ushort e = cast(ushort) ((xe[F.EXPPOS_SHORT] & F.EXPMASK)
+ + (ye[F.EXPPOS_SHORT] & F.EXPMASK));
+ if (m & 0x8000_0000_0000_0000L)
+ {
+ ++e;
+ m &= 0x7FFF_FFFF_FFFF_FFFFL;
+ }
+ // Now do a multi-byte right shift
+ const uint c = e & 1; // carry
+ e >>= 1;
+ m >>>= 1;
+ if (c)
+ m |= 0x4000_0000_0000_0000L; // shift carry into significand
+ if (e)
+ *ul = m | 0x8000_0000_0000_0000L; // set implicit bit...
+ else
+ *ul = m; // ... unless exponent is 0 (subnormal or zero).
+
+ ue[4]= e | (xe[F.EXPPOS_SHORT]& 0x8000); // restore sign bit
+ }
+ else static if (F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ // This would be trivial if 'ucent' were implemented...
+ ulong *ul = cast(ulong *)&u;
+ ulong *xl = cast(ulong *)&x;
+ ulong *yl = cast(ulong *)&y;
+
+ // Multi-byte add, then multi-byte right shift.
+ import core.checkedint : addu;
+ bool carry;
+ ulong ml = addu(xl[MANTISSA_LSB], yl[MANTISSA_LSB], carry);
+
+ ulong mh = carry + (xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL) +
+ (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL);
+
+ ul[MANTISSA_MSB] = (mh >>> 1) | (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000);
+ ul[MANTISSA_LSB] = (ml >>> 1) | (mh & 1) << 63;
+ }
+ else static if (F.realFormat == RealFormat.ieeeDouble)
+ {
+ ulong *ul = cast(ulong *)&u;
+ ulong *xl = cast(ulong *)&x;
+ ulong *yl = cast(ulong *)&y;
+ ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL)
+ + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1;
+ m |= ((*xl) & 0x8000_0000_0000_0000L);
+ *ul = m;
+ }
+ else static if (F.realFormat == RealFormat.ieeeSingle)
+ {
+ uint *ul = cast(uint *)&u;
+ uint *xl = cast(uint *)&x;
+ uint *yl = cast(uint *)&y;
+ uint m = (((*xl) & 0x7FFF_FFFF) + ((*yl) & 0x7FFF_FFFF)) >>> 1;
+ m |= ((*xl) & 0x8000_0000);
+ *ul = m;
+ }
+ else
+ {
+ assert(0, "Not implemented");
+ }
+ return u;
+}
+
+@safe pure nothrow @nogc unittest
+{
+ assert(ieeeMean(-0.0,-1e-20)<0);
+ assert(ieeeMean(0.0,1e-20)>0);
+
+ assert(ieeeMean(1.0L,4.0L)==2L);
+ assert(ieeeMean(2.0*1.013,8.0*1.013)==4*1.013);
+ assert(ieeeMean(-1.0L,-4.0L)==-2L);
+ assert(ieeeMean(-1.0,-4.0)==-2);
+ assert(ieeeMean(-1.0f,-4.0f)==-2f);
+ assert(ieeeMean(-1.0,-2.0)==-1.5);
+ assert(ieeeMean(-1*(1+8*real.epsilon),-2*(1+8*real.epsilon))
+ ==-1.5*(1+5*real.epsilon));
+ assert(ieeeMean(0x1p60,0x1p-10)==0x1p25);
+
+ static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
+ {
+ assert(ieeeMean(1.0L,real.infinity)==0x1p8192L);
+ assert(ieeeMean(0.0L,real.infinity)==1.5);
+ }
+ assert(ieeeMean(0.5*real.min_normal*(1-4*real.epsilon),0.5*real.min_normal)
+ == 0.5*real.min_normal*(1-2*real.epsilon));
+}
+
+public:
+
+
+/***********************************
+ * Evaluate polynomial A(x) = $(SUB a, 0) + $(SUB a, 1)x + $(SUB a, 2)$(POWER x,2)
+ * + $(SUB a,3)$(POWER x,3); ...
+ *
+ * Uses Horner's rule A(x) = $(SUB a, 0) + x($(SUB a, 1) + x($(SUB a, 2)
+ * + x($(SUB a, 3) + ...)))
+ * Params:
+ * x = the value to evaluate.
+ * A = array of coefficients $(SUB a, 0), $(SUB a, 1), etc.
+ */
+Unqual!(CommonType!(T1, T2)) poly(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc
+if (isFloatingPoint!T1 && isFloatingPoint!T2)
+in
+{
+ assert(A.length > 0);
+}
+body
+{
+ static if (is(Unqual!T2 == real))
+ {
+ return polyImpl(x, A);
+ }
+ else
+ {
+ return polyImplBase(x, A);
+ }
+}
+
+///
+@safe nothrow @nogc unittest
+{
+ real x = 3.1;
+ static real[] pp = [56.1, 32.7, 6];
+
+ assert(poly(x, pp) == (56.1L + (32.7L + 6.0L * x) * x));
+}
+
+@safe nothrow @nogc unittest
+{
+ double x = 3.1;
+ static double[] pp = [56.1, 32.7, 6];
+ double y = x;
+ y *= 6.0;
+ y += 32.7;
+ y *= x;
+ y += 56.1;
+ assert(poly(x, pp) == y);
+}
+
+@safe unittest
+{
+ static assert(poly(3.0, [1.0, 2.0, 3.0]) == 34);
+}
+
+private Unqual!(CommonType!(T1, T2)) polyImplBase(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc
+if (isFloatingPoint!T1 && isFloatingPoint!T2)
+{
+ ptrdiff_t i = A.length - 1;
+ typeof(return) r = A[i];
+ while (--i >= 0)
+ {
+ r *= x;
+ r += A[i];
+ }
+ return r;
+}
+
+private real polyImpl(real x, in real[] A) @trusted pure nothrow @nogc
+{
+ version (D_InlineAsm_X86)
+ {
+ if (__ctfe)
+ {
+ return polyImplBase(x, A);
+ }
+ version (Windows)
+ {
+ // BUG: This code assumes a frame pointer in EBP.
+ asm pure nothrow @nogc // assembler by W. Bright
+ {
+ // EDX = (A.length - 1) * real.sizeof
+ mov ECX,A[EBP] ; // ECX = A.length
+ dec ECX ;
+ lea EDX,[ECX][ECX*8] ;
+ add EDX,ECX ;
+ add EDX,A+4[EBP] ;
+ fld real ptr [EDX] ; // ST0 = coeff[ECX]
+ jecxz return_ST ;
+ fld x[EBP] ; // ST0 = x
+ fxch ST(1) ; // ST1 = x, ST0 = r
+ align 4 ;
+ L2: fmul ST,ST(1) ; // r *= x
+ fld real ptr -10[EDX] ;
+ sub EDX,10 ; // deg--
+ faddp ST(1),ST ;
+ dec ECX ;
+ jne L2 ;
+ fxch ST(1) ; // ST1 = r, ST0 = x
+ fstp ST(0) ; // dump x
+ align 4 ;
+ return_ST: ;
+ ;
+ }
+ }
+ else version (linux)
+ {
+ asm pure nothrow @nogc // assembler by W. Bright
+ {
+ // EDX = (A.length - 1) * real.sizeof
+ mov ECX,A[EBP] ; // ECX = A.length
+ dec ECX ;
+ lea EDX,[ECX*8] ;
+ lea EDX,[EDX][ECX*4] ;
+ add EDX,A+4[EBP] ;
+ fld real ptr [EDX] ; // ST0 = coeff[ECX]
+ jecxz return_ST ;
+ fld x[EBP] ; // ST0 = x
+ fxch ST(1) ; // ST1 = x, ST0 = r
+ align 4 ;
+ L2: fmul ST,ST(1) ; // r *= x
+ fld real ptr -12[EDX] ;
+ sub EDX,12 ; // deg--
+ faddp ST(1),ST ;
+ dec ECX ;
+ jne L2 ;
+ fxch ST(1) ; // ST1 = r, ST0 = x
+ fstp ST(0) ; // dump x
+ align 4 ;
+ return_ST: ;
+ ;
+ }
+ }
+ else version (OSX)
+ {
+ asm pure nothrow @nogc // assembler by W. Bright
+ {
+ // EDX = (A.length - 1) * real.sizeof
+ mov ECX,A[EBP] ; // ECX = A.length
+ dec ECX ;
+ lea EDX,[ECX*8] ;
+ add EDX,EDX ;
+ add EDX,A+4[EBP] ;
+ fld real ptr [EDX] ; // ST0 = coeff[ECX]
+ jecxz return_ST ;
+ fld x[EBP] ; // ST0 = x
+ fxch ST(1) ; // ST1 = x, ST0 = r
+ align 4 ;
+ L2: fmul ST,ST(1) ; // r *= x
+ fld real ptr -16[EDX] ;
+ sub EDX,16 ; // deg--
+ faddp ST(1),ST ;
+ dec ECX ;
+ jne L2 ;
+ fxch ST(1) ; // ST1 = r, ST0 = x
+ fstp ST(0) ; // dump x
+ align 4 ;
+ return_ST: ;
+ ;
+ }
+ }
+ else version (FreeBSD)
+ {
+ asm pure nothrow @nogc // assembler by W. Bright
+ {
+ // EDX = (A.length - 1) * real.sizeof
+ mov ECX,A[EBP] ; // ECX = A.length
+ dec ECX ;
+ lea EDX,[ECX*8] ;
+ lea EDX,[EDX][ECX*4] ;
+ add EDX,A+4[EBP] ;
+ fld real ptr [EDX] ; // ST0 = coeff[ECX]
+ jecxz return_ST ;
+ fld x[EBP] ; // ST0 = x
+ fxch ST(1) ; // ST1 = x, ST0 = r
+ align 4 ;
+ L2: fmul ST,ST(1) ; // r *= x
+ fld real ptr -12[EDX] ;
+ sub EDX,12 ; // deg--
+ faddp ST(1),ST ;
+ dec ECX ;
+ jne L2 ;
+ fxch ST(1) ; // ST1 = r, ST0 = x
+ fstp ST(0) ; // dump x
+ align 4 ;
+ return_ST: ;
+ ;
+ }
+ }
+ else version (Solaris)
+ {
+ asm pure nothrow @nogc // assembler by W. Bright
+ {
+ // EDX = (A.length - 1) * real.sizeof
+ mov ECX,A[EBP] ; // ECX = A.length
+ dec ECX ;
+ lea EDX,[ECX*8] ;
+ lea EDX,[EDX][ECX*4] ;
+ add EDX,A+4[EBP] ;
+ fld real ptr [EDX] ; // ST0 = coeff[ECX]
+ jecxz return_ST ;
+ fld x[EBP] ; // ST0 = x
+ fxch ST(1) ; // ST1 = x, ST0 = r
+ align 4 ;
+ L2: fmul ST,ST(1) ; // r *= x
+ fld real ptr -12[EDX] ;
+ sub EDX,12 ; // deg--
+ faddp ST(1),ST ;
+ dec ECX ;
+ jne L2 ;
+ fxch ST(1) ; // ST1 = r, ST0 = x
+ fstp ST(0) ; // dump x
+ align 4 ;
+ return_ST: ;
+ ;
+ }
+ }
+ else
+ {
+ static assert(0);
+ }
+ }
+ else
+ {
+ return polyImplBase(x, A);
+ }
+}
+
+
+/**
+ Computes whether two values are approximately equal, admitting a maximum
+ relative difference, and a maximum absolute difference.
+
+ Params:
+ lhs = First item to compare.
+ rhs = Second item to compare.
+ maxRelDiff = Maximum allowable difference relative to `rhs`.
+ maxAbsDiff = Maximum absolute difference.
+
+ Returns:
+ `true` if the two items are approximately equal under either criterium.
+ If one item is a range, and the other is a single value, then the result
+ is the logical and-ing of calling `approxEqual` on each element of the
+ ranged item against the single item. If both items are ranges, then
+ `approxEqual` returns `true` if and only if the ranges have the same
+ number of elements and if `approxEqual` evaluates to `true` for each
+ pair of elements.
+ */
+bool approxEqual(T, U, V)(T lhs, U rhs, V maxRelDiff, V maxAbsDiff = 1e-5)
+{
+ import std.range.primitives : empty, front, isInputRange, popFront;
+ static if (isInputRange!T)
+ {
+ static if (isInputRange!U)
+ {
+ // Two ranges
+ for (;; lhs.popFront(), rhs.popFront())
+ {
+ if (lhs.empty) return rhs.empty;
+ if (rhs.empty) return lhs.empty;
+ if (!approxEqual(lhs.front, rhs.front, maxRelDiff, maxAbsDiff))
+ return false;
+ }
+ }
+ else static if (isIntegral!U)
+ {
+ // convert rhs to real
+ return approxEqual(lhs, real(rhs), maxRelDiff, maxAbsDiff);
+ }
+ else
+ {
+ // lhs is range, rhs is number
+ for (; !lhs.empty; lhs.popFront())
+ {
+ if (!approxEqual(lhs.front, rhs, maxRelDiff, maxAbsDiff))
+ return false;
+ }
+ return true;
+ }
+ }
+ else
+ {
+ static if (isInputRange!U)
+ {
+ // lhs is number, rhs is range
+ for (; !rhs.empty; rhs.popFront())
+ {
+ if (!approxEqual(lhs, rhs.front, maxRelDiff, maxAbsDiff))
+ return false;
+ }
+ return true;
+ }
+ else static if (isIntegral!T || isIntegral!U)
+ {
+ // convert both lhs and rhs to real
+ return approxEqual(real(lhs), real(rhs), maxRelDiff, maxAbsDiff);
+ }
+ else
+ {
+ // two numbers
+ //static assert(is(T : real) && is(U : real));
+ if (rhs == 0)
+ {
+ return fabs(lhs) <= maxAbsDiff;
+ }
+ static if (is(typeof(lhs.infinity)) && is(typeof(rhs.infinity)))
+ {
+ if (lhs == lhs.infinity && rhs == rhs.infinity ||
+ lhs == -lhs.infinity && rhs == -rhs.infinity) return true;
+ }
+ return fabs((lhs - rhs) / rhs) <= maxRelDiff
+ || maxAbsDiff != 0 && fabs(lhs - rhs) <= maxAbsDiff;
+ }
+ }
+}
+
+/**
+ Returns $(D approxEqual(lhs, rhs, 1e-2, 1e-5)).
+ */
+bool approxEqual(T, U)(T lhs, U rhs)
+{
+ return approxEqual(lhs, rhs, 1e-2, 1e-5);
+}
+
+///
+@safe pure nothrow unittest
+{
+ assert(approxEqual(1.0, 1.0099));
+ assert(!approxEqual(1.0, 1.011));
+ float[] arr1 = [ 1.0, 2.0, 3.0 ];
+ double[] arr2 = [ 1.001, 1.999, 3 ];
+ assert(approxEqual(arr1, arr2));
+
+ real num = real.infinity;
+ assert(num == real.infinity); // Passes.
+ assert(approxEqual(num, real.infinity)); // Fails.
+ num = -real.infinity;
+ assert(num == -real.infinity); // Passes.
+ assert(approxEqual(num, -real.infinity)); // Fails.
+
+ assert(!approxEqual(3, 0));
+ assert(approxEqual(3, 3));
+ assert(approxEqual(3.0, 3));
+ assert(approxEqual([3, 3, 3], 3.0));
+ assert(approxEqual([3.0, 3.0, 3.0], 3));
+ int a = 10;
+ assert(approxEqual(10, a));
+}
+
+@safe pure nothrow @nogc unittest
+{
+ real num = real.infinity;
+ assert(num == real.infinity); // Passes.
+ assert(approxEqual(num, real.infinity)); // Fails.
+}
+
+
+@safe pure nothrow @nogc unittest
+{
+ float f = sqrt(2.0f);
+ assert(fabs(f * f - 2.0f) < .00001);
+
+ double d = sqrt(2.0);
+ assert(fabs(d * d - 2.0) < .00001);
+
+ real r = sqrt(2.0L);
+ assert(fabs(r * r - 2.0) < .00001);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ float f = fabs(-2.0f);
+ assert(f == 2);
+
+ double d = fabs(-2.0);
+ assert(d == 2);
+
+ real r = fabs(-2.0L);
+ assert(r == 2);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ float f = sin(-2.0f);
+ assert(fabs(f - -0.909297f) < .00001);
+
+ double d = sin(-2.0);
+ assert(fabs(d - -0.909297f) < .00001);
+
+ real r = sin(-2.0L);
+ assert(fabs(r - -0.909297f) < .00001);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ float f = cos(-2.0f);
+ assert(fabs(f - -0.416147f) < .00001);
+
+ double d = cos(-2.0);
+ assert(fabs(d - -0.416147f) < .00001);
+
+ real r = cos(-2.0L);
+ assert(fabs(r - -0.416147f) < .00001);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ float f = tan(-2.0f);
+ assert(fabs(f - 2.18504f) < .00001);
+
+ double d = tan(-2.0);
+ assert(fabs(d - 2.18504f) < .00001);
+
+ real r = tan(-2.0L);
+ assert(fabs(r - 2.18504f) < .00001);
+
+ // Verify correct behavior for large inputs
+ assert(!isNaN(tan(0x1p63)));
+ assert(!isNaN(tan(0x1p300L)));
+ assert(!isNaN(tan(-0x1p63)));
+ assert(!isNaN(tan(-0x1p300L)));
+}
+
+@safe pure nothrow unittest
+{
+ // issue 6381: floor/ceil should be usable in pure function.
+ auto x = floor(1.2);
+ auto y = ceil(1.2);
+}
+
+@safe pure nothrow unittest
+{
+ // relative comparison depends on rhs, make sure proper side is used when
+ // comparing range to single value. Based on bugzilla issue 15763
+ auto a = [2e-3 - 1e-5];
+ auto b = 2e-3 + 1e-5;
+ assert(a[0].approxEqual(b));
+ assert(!b.approxEqual(a[0]));
+ assert(a.approxEqual(b));
+ assert(!b.approxEqual(a));
+}
+
+/***********************************
+ * Defines a total order on all floating-point numbers.
+ *
+ * The order is defined as follows:
+ * $(UL
+ * $(LI All numbers in [-$(INFIN), +$(INFIN)] are ordered
+ * the same way as by built-in comparison, with the exception of
+ * -0.0, which is less than +0.0;)
+ * $(LI If the sign bit is set (that is, it's 'negative'), $(NAN) is less
+ * than any number; if the sign bit is not set (it is 'positive'),
+ * $(NAN) is greater than any number;)
+ * $(LI $(NAN)s of the same sign are ordered by the payload ('negative'
+ * ones - in reverse order).)
+ * )
+ *
+ * Returns:
+ * negative value if $(D x) precedes $(D y) in the order specified above;
+ * 0 if $(D x) and $(D y) are identical, and positive value otherwise.
+ *
+ * See_Also:
+ * $(MYREF isIdentical)
+ * Standards: Conforms to IEEE 754-2008
+ */
+int cmp(T)(const(T) x, const(T) y) @nogc @trusted pure nothrow
+if (isFloatingPoint!T)
+{
+ alias F = floatTraits!T;
+
+ static if (F.realFormat == RealFormat.ieeeSingle
+ || F.realFormat == RealFormat.ieeeDouble)
+ {
+ static if (T.sizeof == 4)
+ alias UInt = uint;
+ else
+ alias UInt = ulong;
+
+ union Repainter
+ {
+ T number;
+ UInt bits;
+ }
+
+ enum msb = ~(UInt.max >>> 1);
+
+ import std.typecons : Tuple;
+ Tuple!(Repainter, Repainter) vars = void;
+ vars[0].number = x;
+ vars[1].number = y;
+
+ foreach (ref var; vars)
+ if (var.bits & msb)
+ var.bits = ~var.bits;
+ else
+ var.bits |= msb;
+
+ if (vars[0].bits < vars[1].bits)
+ return -1;
+ else if (vars[0].bits > vars[1].bits)
+ return 1;
+ else
+ return 0;
+ }
+ else static if (F.realFormat == RealFormat.ieeeExtended53
+ || F.realFormat == RealFormat.ieeeExtended
+ || F.realFormat == RealFormat.ieeeQuadruple)
+ {
+ static if (F.realFormat == RealFormat.ieeeQuadruple)
+ alias RemT = ulong;
+ else
+ alias RemT = ushort;
+
+ struct Bits
+ {
+ ulong bulk;
+ RemT rem;
+ }
+
+ union Repainter
+ {
+ T number;
+ Bits bits;
+ ubyte[T.sizeof] bytes;
+ }
+
+ import std.typecons : Tuple;
+ Tuple!(Repainter, Repainter) vars = void;
+ vars[0].number = x;
+ vars[1].number = y;
+
+ foreach (ref var; vars)
+ if (var.bytes[F.SIGNPOS_BYTE] & 0x80)
+ {
+ var.bits.bulk = ~var.bits.bulk;
+ var.bits.rem = cast(typeof(var.bits.rem))(-1 - var.bits.rem); // ~var.bits.rem
+ }
+ else
+ {
+ var.bytes[F.SIGNPOS_BYTE] |= 0x80;
+ }
+
+ version (LittleEndian)
+ {
+ if (vars[0].bits.rem < vars[1].bits.rem)
+ return -1;
+ else if (vars[0].bits.rem > vars[1].bits.rem)
+ return 1;
+ else if (vars[0].bits.bulk < vars[1].bits.bulk)
+ return -1;
+ else if (vars[0].bits.bulk > vars[1].bits.bulk)
+ return 1;
+ else
+ return 0;
+ }
+ else
+ {
+ if (vars[0].bits.bulk < vars[1].bits.bulk)
+ return -1;
+ else if (vars[0].bits.bulk > vars[1].bits.bulk)
+ return 1;
+ else if (vars[0].bits.rem < vars[1].bits.rem)
+ return -1;
+ else if (vars[0].bits.rem > vars[1].bits.rem)
+ return 1;
+ else
+ return 0;
+ }
+ }
+ else
+ {
+ // IBM Extended doubledouble does not follow the general
+ // sign-exponent-significand layout, so has to be handled generically
+
+ const int xSign = signbit(x),
+ ySign = signbit(y);
+
+ if (xSign == 1 && ySign == 1)
+ return cmp(-y, -x);
+ else if (xSign == 1)
+ return -1;
+ else if (ySign == 1)
+ return 1;
+ else if (x < y)
+ return -1;
+ else if (x == y)
+ return 0;
+ else if (x > y)
+ return 1;
+ else if (isNaN(x) && !isNaN(y))
+ return 1;
+ else if (isNaN(y) && !isNaN(x))
+ return -1;
+ else if (getNaNPayload(x) < getNaNPayload(y))
+ return -1;
+ else if (getNaNPayload(x) > getNaNPayload(y))
+ return 1;
+ else
+ return 0;
+ }
+}
+
+/// Most numbers are ordered naturally.
+@safe unittest
+{
+ assert(cmp(-double.infinity, -double.max) < 0);
+ assert(cmp(-double.max, -100.0) < 0);
+ assert(cmp(-100.0, -0.5) < 0);
+ assert(cmp(-0.5, 0.0) < 0);
+ assert(cmp(0.0, 0.5) < 0);
+ assert(cmp(0.5, 100.0) < 0);
+ assert(cmp(100.0, double.max) < 0);
+ assert(cmp(double.max, double.infinity) < 0);
+
+ assert(cmp(1.0, 1.0) == 0);
+}
+
+/// Positive and negative zeroes are distinct.
+@safe unittest
+{
+ assert(cmp(-0.0, +0.0) < 0);
+ assert(cmp(+0.0, -0.0) > 0);
+}
+
+/// Depending on the sign, $(NAN)s go to either end of the spectrum.
+@safe unittest
+{
+ assert(cmp(-double.nan, -double.infinity) < 0);
+ assert(cmp(double.infinity, double.nan) < 0);
+ assert(cmp(-double.nan, double.nan) < 0);
+}
+
+/// $(NAN)s of the same sign are ordered by the payload.
+@safe unittest
+{
+ assert(cmp(NaN(10), NaN(20)) < 0);
+ assert(cmp(-NaN(20), -NaN(10)) < 0);
+}
+
+@safe unittest
+{
+ import std.meta : AliasSeq;
+ foreach (T; AliasSeq!(float, double, real))
+ {
+ T[] values = [-cast(T) NaN(20), -cast(T) NaN(10), -T.nan, -T.infinity,
+ -T.max, -T.max / 2, T(-16.0), T(-1.0).nextDown,
+ T(-1.0), T(-1.0).nextUp,
+ T(-0.5), -T.min_normal, (-T.min_normal).nextUp,
+ -2 * T.min_normal * T.epsilon,
+ -T.min_normal * T.epsilon,
+ T(-0.0), T(0.0),
+ T.min_normal * T.epsilon,
+ 2 * T.min_normal * T.epsilon,
+ T.min_normal.nextDown, T.min_normal, T(0.5),
+ T(1.0).nextDown, T(1.0),
+ T(1.0).nextUp, T(16.0), T.max / 2, T.max,
+ T.infinity, T.nan, cast(T) NaN(10), cast(T) NaN(20)];
+
+ foreach (i, x; values)
+ {
+ foreach (y; values[i + 1 .. $])
+ {
+ assert(cmp(x, y) < 0);
+ assert(cmp(y, x) > 0);
+ }
+ assert(cmp(x, x) == 0);
+ }
+ }
+}
+
+private enum PowType
+{
+ floor,
+ ceil
+}
+
+pragma(inline, true)
+private T powIntegralImpl(PowType type, T)(T val)
+{
+ import core.bitop : bsr;
+
+ if (val == 0 || (type == PowType.ceil && (val > T.max / 2 || val == T.min)))
+ return 0;
+ else
+ {
+ static if (isSigned!T)
+ return cast(Unqual!T) (val < 0 ? -(T(1) << bsr(0 - val) + type) : T(1) << bsr(val) + type);
+ else
+ return cast(Unqual!T) (T(1) << bsr(val) + type);
+ }
+}
+
+private T powFloatingPointImpl(PowType type, T)(T x)
+{
+ if (!x.isFinite)
+ return x;
+
+ if (!x)
+ return x;
+
+ int exp;
+ auto y = frexp(x, exp);
+
+ static if (type == PowType.ceil)
+ y = ldexp(cast(T) 0.5, exp + 1);
+ else
+ y = ldexp(cast(T) 0.5, exp);
+
+ if (!y.isFinite)
+ return cast(T) 0.0;
+
+ y = copysign(y, x);
+
+ return y;
+}
+
+/**
+ * Gives the next power of two after $(D val). `T` can be any built-in
+ * numerical type.
+ *
+ * If the operation would lead to an over/underflow, this function will
+ * return `0`.
+ *
+ * Params:
+ * val = any number
+ *
+ * Returns:
+ * the next power of two after $(D val)
+ */
+T nextPow2(T)(const T val)
+if (isIntegral!T)
+{
+ return powIntegralImpl!(PowType.ceil)(val);
+}
+
+/// ditto
+T nextPow2(T)(const T val)
+if (isFloatingPoint!T)
+{
+ return powFloatingPointImpl!(PowType.ceil)(val);
+}
+
+///
+@safe @nogc pure nothrow unittest
+{
+ assert(nextPow2(2) == 4);
+ assert(nextPow2(10) == 16);
+ assert(nextPow2(4000) == 4096);
+
+ assert(nextPow2(-2) == -4);
+ assert(nextPow2(-10) == -16);
+
+ assert(nextPow2(uint.max) == 0);
+ assert(nextPow2(uint.min) == 0);
+ assert(nextPow2(size_t.max) == 0);
+ assert(nextPow2(size_t.min) == 0);
+
+ assert(nextPow2(int.max) == 0);
+ assert(nextPow2(int.min) == 0);
+ assert(nextPow2(long.max) == 0);
+ assert(nextPow2(long.min) == 0);
+}
+
+///
+@safe @nogc pure nothrow unittest
+{
+ assert(nextPow2(2.1) == 4.0);
+ assert(nextPow2(-2.0) == -4.0);
+ assert(nextPow2(0.25) == 0.5);
+ assert(nextPow2(-4.0) == -8.0);
+
+ assert(nextPow2(double.max) == 0.0);
+ assert(nextPow2(double.infinity) == double.infinity);
+}
+
+@safe @nogc pure nothrow unittest
+{
+ assert(nextPow2(ubyte(2)) == 4);
+ assert(nextPow2(ubyte(10)) == 16);
+
+ assert(nextPow2(byte(2)) == 4);
+ assert(nextPow2(byte(10)) == 16);
+
+ assert(nextPow2(short(2)) == 4);
+ assert(nextPow2(short(10)) == 16);
+ assert(nextPow2(short(4000)) == 4096);
+
+ assert(nextPow2(ushort(2)) == 4);
+ assert(nextPow2(ushort(10)) == 16);
+ assert(nextPow2(ushort(4000)) == 4096);
+}
+
+@safe @nogc pure nothrow unittest
+{
+ foreach (ulong i; 1 .. 62)
+ {
+ assert(nextPow2(1UL << i) == 2UL << i);
+ assert(nextPow2((1UL << i) - 1) == 1UL << i);
+ assert(nextPow2((1UL << i) + 1) == 2UL << i);
+ assert(nextPow2((1UL << i) + (1UL<<(i-1))) == 2UL << i);
+ }
+}
+
+@safe @nogc pure nothrow unittest
+{
+ import std.meta : AliasSeq;
+
+ foreach (T; AliasSeq!(float, double, real))
+ {
+ enum T subNormal = T.min_normal / 2;
+
+ static if (subNormal) assert(nextPow2(subNormal) == T.min_normal);
+
+ assert(nextPow2(T(0.0)) == 0.0);
+
+ assert(nextPow2(T(2.0)) == 4.0);
+ assert(nextPow2(T(2.1)) == 4.0);
+ assert(nextPow2(T(3.1)) == 4.0);
+ assert(nextPow2(T(4.0)) == 8.0);
+ assert(nextPow2(T(0.25)) == 0.5);
+
+ assert(nextPow2(T(-2.0)) == -4.0);
+ assert(nextPow2(T(-2.1)) == -4.0);
+ assert(nextPow2(T(-3.1)) == -4.0);
+ assert(nextPow2(T(-4.0)) == -8.0);
+ assert(nextPow2(T(-0.25)) == -0.5);
+
+ assert(nextPow2(T.max) == 0);
+ assert(nextPow2(-T.max) == 0);
+
+ assert(nextPow2(T.infinity) == T.infinity);
+ assert(nextPow2(T.init).isNaN);
+ }
+}
+
+@safe @nogc pure nothrow unittest // Issue 15973
+{
+ assert(nextPow2(uint.max / 2) == uint.max / 2 + 1);
+ assert(nextPow2(uint.max / 2 + 2) == 0);
+ assert(nextPow2(int.max / 2) == int.max / 2 + 1);
+ assert(nextPow2(int.max / 2 + 2) == 0);
+ assert(nextPow2(int.min + 1) == int.min);
+}
+
+/**
+ * Gives the last power of two before $(D val). $(T) can be any built-in
+ * numerical type.
+ *
+ * Params:
+ * val = any number
+ *
+ * Returns:
+ * the last power of two before $(D val)
+ */
+T truncPow2(T)(const T val)
+if (isIntegral!T)
+{
+ return powIntegralImpl!(PowType.floor)(val);
+}
+
+/// ditto
+T truncPow2(T)(const T val)
+if (isFloatingPoint!T)
+{
+ return powFloatingPointImpl!(PowType.floor)(val);
+}
+
+///
+@safe @nogc pure nothrow unittest
+{
+ assert(truncPow2(3) == 2);
+ assert(truncPow2(4) == 4);
+ assert(truncPow2(10) == 8);
+ assert(truncPow2(4000) == 2048);
+
+ assert(truncPow2(-5) == -4);
+ assert(truncPow2(-20) == -16);
+
+ assert(truncPow2(uint.max) == int.max + 1);
+ assert(truncPow2(uint.min) == 0);
+ assert(truncPow2(ulong.max) == long.max + 1);
+ assert(truncPow2(ulong.min) == 0);
+
+ assert(truncPow2(int.max) == (int.max / 2) + 1);
+ assert(truncPow2(int.min) == int.min);
+ assert(truncPow2(long.max) == (long.max / 2) + 1);
+ assert(truncPow2(long.min) == long.min);
+}
+
+///
+@safe @nogc pure nothrow unittest
+{
+ assert(truncPow2(2.1) == 2.0);
+ assert(truncPow2(7.0) == 4.0);
+ assert(truncPow2(-1.9) == -1.0);
+ assert(truncPow2(0.24) == 0.125);
+ assert(truncPow2(-7.0) == -4.0);
+
+ assert(truncPow2(double.infinity) == double.infinity);
+}
+
+@safe @nogc pure nothrow unittest
+{
+ assert(truncPow2(ubyte(3)) == 2);
+ assert(truncPow2(ubyte(4)) == 4);
+ assert(truncPow2(ubyte(10)) == 8);
+
+ assert(truncPow2(byte(3)) == 2);
+ assert(truncPow2(byte(4)) == 4);
+ assert(truncPow2(byte(10)) == 8);
+
+ assert(truncPow2(ushort(3)) == 2);
+ assert(truncPow2(ushort(4)) == 4);
+ assert(truncPow2(ushort(10)) == 8);
+ assert(truncPow2(ushort(4000)) == 2048);
+
+ assert(truncPow2(short(3)) == 2);
+ assert(truncPow2(short(4)) == 4);
+ assert(truncPow2(short(10)) == 8);
+ assert(truncPow2(short(4000)) == 2048);
+}
+
+@safe @nogc pure nothrow unittest
+{
+ foreach (ulong i; 1 .. 62)
+ {
+ assert(truncPow2(2UL << i) == 2UL << i);
+ assert(truncPow2((2UL << i) + 1) == 2UL << i);
+ assert(truncPow2((2UL << i) - 1) == 1UL << i);
+ assert(truncPow2((2UL << i) - (2UL<<(i-1))) == 1UL << i);
+ }
+}
+
+@safe @nogc pure nothrow unittest
+{
+ import std.meta : AliasSeq;
+
+ foreach (T; AliasSeq!(float, double, real))
+ {
+ assert(truncPow2(T(0.0)) == 0.0);
+
+ assert(truncPow2(T(4.0)) == 4.0);
+ assert(truncPow2(T(2.1)) == 2.0);
+ assert(truncPow2(T(3.5)) == 2.0);
+ assert(truncPow2(T(7.0)) == 4.0);
+ assert(truncPow2(T(0.24)) == 0.125);
+
+ assert(truncPow2(T(-2.0)) == -2.0);
+ assert(truncPow2(T(-2.1)) == -2.0);
+ assert(truncPow2(T(-3.1)) == -2.0);
+ assert(truncPow2(T(-7.0)) == -4.0);
+ assert(truncPow2(T(-0.24)) == -0.125);
+
+ assert(truncPow2(T.infinity) == T.infinity);
+ assert(truncPow2(T.init).isNaN);
+ }
+}
+
+/**
+Check whether a number is an integer power of two.
+
+Note that only positive numbers can be integer powers of two. This
+function always return `false` if `x` is negative or zero.
+
+Params:
+ x = the number to test
+
+Returns:
+ `true` if `x` is an integer power of two.
+*/
+bool isPowerOf2(X)(const X x) pure @safe nothrow @nogc
+if (isNumeric!X)
+{
+ static if (isFloatingPoint!X)
+ {
+ int exp;
+ const X sig = frexp(x, exp);
+
+ return (exp != int.min) && (sig is cast(X) 0.5L);
+ }
+ else
+ {
+ static if (isSigned!X)
+ {
+ auto y = cast(typeof(x + 0))x;
+ return y > 0 && !(y & (y - 1));
+ }
+ else
+ {
+ auto y = cast(typeof(x + 0u))x;
+ return (y & -y) > (y - 1);
+ }
+ }
+}
+///
+@safe unittest
+{
+ assert( isPowerOf2(1.0L));
+ assert( isPowerOf2(2.0L));
+ assert( isPowerOf2(0.5L));
+ assert( isPowerOf2(pow(2.0L, 96)));
+ assert( isPowerOf2(pow(2.0L, -77)));
+
+ assert(!isPowerOf2(-2.0L));
+ assert(!isPowerOf2(-0.5L));
+ assert(!isPowerOf2(0.0L));
+ assert(!isPowerOf2(4.315));
+ assert(!isPowerOf2(1.0L / 3.0L));
+
+ assert(!isPowerOf2(real.nan));
+ assert(!isPowerOf2(real.infinity));
+}
+///
+@safe unittest
+{
+ assert( isPowerOf2(1));
+ assert( isPowerOf2(2));
+ assert( isPowerOf2(1uL << 63));
+
+ assert(!isPowerOf2(-4));
+ assert(!isPowerOf2(0));
+ assert(!isPowerOf2(1337u));
+}
+
+@safe unittest
+{
+ import std.meta : AliasSeq;
+
+ immutable smallP2 = pow(2.0L, -62);
+ immutable bigP2 = pow(2.0L, 50);
+ immutable smallP7 = pow(7.0L, -35);
+ immutable bigP7 = pow(7.0L, 30);
+
+ foreach (X; AliasSeq!(float, double, real))
+ {
+ immutable min_sub = X.min_normal * X.epsilon;
+
+ foreach (x; AliasSeq!(smallP2, min_sub, X.min_normal, .25L, 0.5L, 1.0L,
+ 2.0L, 8.0L, pow(2.0L, X.max_exp - 1), bigP2))
+ {
+ assert( isPowerOf2(cast(X) x));
+ assert(!isPowerOf2(cast(X)-x));
+ }
+
+ foreach (x; AliasSeq!(0.0L, 3 * min_sub, smallP7, 0.1L, 1337.0L, bigP7, X.max, real.nan, real.infinity))
+ {
+ assert(!isPowerOf2(cast(X) x));
+ assert(!isPowerOf2(cast(X)-x));
+ }
+ }
+
+ foreach (X; AliasSeq!(byte, ubyte, short, ushort, int, uint, long, ulong))
+ {
+ foreach (x; [1, 2, 4, 8, (X.max >>> 1) + 1])
+ {
+ assert( isPowerOf2(cast(X) x));
+ static if (isSigned!X)
+ assert(!isPowerOf2(cast(X)-x));
+ }
+
+ foreach (x; [0, 3, 5, 13, 77, X.min, X.max])
+ assert(!isPowerOf2(cast(X) x));
+ }
+}
generated by cgit v1.1 at 2025-03-31 15:59:20 +0000