summaryrefslogtreecommitdiff
path: root/libphobos/src/std/math/rounding.d
diff options
context:
space:
mode:
Diffstat (limited to 'libphobos/src/std/math/rounding.d')
-rw-r--r--libphobos/src/std/math/rounding.d1004
1 files changed, 1004 insertions, 0 deletions
diff --git a/libphobos/src/std/math/rounding.d b/libphobos/src/std/math/rounding.d
new file mode 100644
index 00000000000..5c8d708c489
--- /dev/null
+++ b/libphobos/src/std/math/rounding.d
@@ -0,0 +1,1004 @@
+// Written in the D programming language.
+
+/**
+This is a submodule of $(MREF std, math).
+
+It contains several functions for rounding floating point numbers.
+
+Copyright: Copyright The D Language Foundation 2000 - 2011.
+ D implementations of 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/rounding.d)
+ */
+
+/* NOTE: This file has been patched from the original DMD distribution to
+ * work with the GDC compiler.
+ */
+module std.math.rounding;
+
+static import core.math;
+static import core.stdc.math;
+
+import std.traits : isFloatingPoint, isIntegral, Unqual;
+
+version (D_InlineAsm_X86) version = InlineAsm_X86_Any;
+version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any;
+
+version (InlineAsm_X86_Any) version = InlineAsm_X87;
+version (InlineAsm_X87)
+{
+ static assert(real.mant_dig == 64);
+ version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC;
+}
+
+/**************************************
+ * Returns the value of x rounded upward to the next integer
+ * (toward positive infinity).
+ */
+real ceil(real x) @trusted pure nothrow @nogc
+{
+ version (InlineAsm_X87_MSVC)
+ {
+ version (X86_64)
+ {
+ 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
+ {
+ 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
+ {
+ import std.math.traits : isInfinity, isNaN;
+
+ // 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
+{
+ import std.math.traits : isNaN;
+
+ 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
+{
+ import std.math.traits : isInfinity, isNaN;
+
+ // 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
+{
+ import std.math.traits : isNaN;
+
+ 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
+{
+ import std.math.traits : isInfinity, isNaN;
+
+ // 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
+{
+ import std.math.traits : isNaN;
+
+ 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 (InlineAsm_X87_MSVC)
+ {
+ version (X86_64)
+ {
+ 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
+ {
+ 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
+ {
+ import std.math.traits : isInfinity, isNaN;
+
+ // Special cases.
+ if (isNaN(x) || isInfinity(x) || x == 0.0)
+ return x;
+
+ return floorImpl(x);
+ }
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ import std.math.traits : isNaN;
+
+ assert(floor(+123.456L) == +123);
+ assert(floor(-123.456L) == -124);
+ assert(floor(+123.0L) == +123);
+ assert(floor(-124.0L) == -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
+{
+ import std.math.traits : isInfinity, isNaN;
+
+ // Special cases.
+ if (isNaN(x) || isInfinity(x) || x == 0.0)
+ return x;
+
+ return floorImpl(x);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ import std.math.traits : isNaN;
+
+ assert(floor(+123.456) == +123);
+ assert(floor(-123.456) == -124);
+ assert(floor(+123.0) == +123);
+ assert(floor(-124.0) == -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
+{
+ import std.math.traits : isInfinity, isNaN;
+
+ // Special cases.
+ if (isNaN(x) || isInfinity(x) || x == 0.0)
+ return x;
+
+ return floorImpl(x);
+}
+
+@safe pure nothrow @nogc unittest
+{
+ import std.math.traits : isNaN;
+
+ assert(floor(+123.456f) == +123);
+ assert(floor(-123.456f) == -124);
+ assert(floor(+123.0f) == +123);
+ assert(floor(-124.0f) == -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)));
+}
+
+// https://issues.dlang.org/show_bug.cgi?id=6381
+// floor/ceil should be usable in pure function.
+@safe pure nothrow unittest
+{
+ auto x = floor(1.2);
+ auto y = ceil(1.2);
+}
+
+/**
+ * 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)
+{
+ import std.math.traits : isInfinity;
+
+ 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
+{
+ import std.math.operations : isClose;
+
+ assert(isClose(12345.6789L.quantize(0.01L), 12345.68L));
+ assert(isClose(12345.6789L.quantize!floor(0.01L), 12345.67L));
+ assert(isClose(12345.6789L.quantize(22.0L), 12342.0L));
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ import std.math.operations : isClose;
+ import std.math.traits : isNaN;
+
+ assert(isClose(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)
+{
+ import std.math.exponential : pow;
+
+ // 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)
+{
+ import std.math.exponential : pow;
+
+ enum unit = cast(F) pow(base, exp);
+ return quantize!rfunc(val, unit);
+}
+
+///
+@safe pure nothrow @nogc unittest
+{
+ import std.math.operations : isClose;
+
+ assert(isClose(12345.6789L.quantize!10(-2), 12345.68L));
+ assert(isClose(12345.6789L.quantize!(10, -2), 12345.68L));
+ assert(isClose(12345.6789L.quantize!(10, floor)(-2), 12345.67L));
+ assert(isClose(12345.6789L.quantize!(10, -2, floor), 12345.67L));
+
+ assert(isClose(12345.6789L.quantize!22(1), 12342.0L));
+ assert(isClose(12345.6789L.quantize!22, 12342.0L));
+}
+
+@safe pure nothrow @nogc unittest
+{
+ import std.math.exponential : log10, pow;
+ import std.math.operations : isClose;
+ import std.meta : AliasSeq;
+
+ static foreach (F; AliasSeq!(real, double, float))
+ {{
+ const maxL10 = cast(int) F.max.log10.floor;
+ const maxR10 = pow(cast(F) 10, maxL10);
+ assert(isClose((cast(F) 0.9L * maxR10).quantize!10(maxL10), maxR10));
+ assert(isClose((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.
+ */
+pragma(inline, true)
+real nearbyint(real x) @safe pure nothrow @nogc
+{
+ return core.stdc.math.nearbyintl(x);
+}
+
+///
+@safe pure unittest
+{
+ import std.math.traits : isNaN;
+
+ assert(nearbyint(0.4) == 0);
+ assert(nearbyint(0.5) == 0);
+ assert(nearbyint(0.6) == 1);
+ assert(nearbyint(100.0) == 100);
+
+ assert(isNaN(nearbyint(real.nan)));
+ assert(nearbyint(real.infinity) == real.infinity);
+ assert(nearbyint(-real.infinity) == -real.infinity);
+}
+
+/**********************************
+ * 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.
+ *
+ * $(LREF nearbyint) performs the same operation, but does
+ * not set the FE_INEXACT exception.
+ */
+pragma(inline, true)
+real rint(real x) @safe pure nothrow @nogc
+{
+ return core.math.rint(x);
+}
+///ditto
+pragma(inline, true)
+double rint(double x) @safe pure nothrow @nogc
+{
+ return core.math.rint(x);
+}
+///ditto
+pragma(inline, true)
+float rint(float x) @safe pure nothrow @nogc
+{
+ return core.math.rint(x);
+}
+
+///
+@safe unittest
+{
+ import std.math.traits : isNaN;
+
+ version (IeeeFlagsSupport) resetIeeeFlags();
+ assert(rint(0.4) == 0);
+ version (GNU) { /* inexact bit not set with enabled optimizations */ } else
+ version (IeeeFlagsSupport) assert(ieeeFlags.inexact);
+
+ assert(rint(0.5) == 0);
+ assert(rint(0.6) == 1);
+ assert(rint(100.0) == 100);
+
+ assert(isNaN(rint(real.nan)));
+ assert(rint(real.infinity) == real.infinity);
+ assert(rint(-real.infinity) == -real.infinity);
+}
+
+@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_X87)
+ {
+ 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
+ {
+ import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
+
+ 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 ||
+ F.realFormat == RealFormat.ieeeExtended53)
+ {
+ long result;
+
+ // Rounding limit when casting from real(80-bit) to ulong.
+ static if (F.realFormat == RealFormat.ieeeExtended)
+ enum real OF = 9.22337203685477580800E18L;
+ else
+ enum real OF = 4.50359962737049600000E15L;
+
+ 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 exp = (vu[F.EXPPOS_SHORT] & F.EXPMASK) - (F.EXPBIAS + 1);
+ const implicitOne = 1UL << 48;
+ auto vl = cast(ulong*)(&x);
+ vl[MANTISSA_MSB] &= implicitOne - 1;
+ vl[MANTISSA_MSB] |= implicitOne;
+
+ long result;
+
+ 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.
+ *
+ * Returns:
+ * A `real`.
+ */
+auto round(real x) @trusted nothrow @nogc
+{
+ version (CRuntime_Microsoft)
+ {
+ import std.math.hardware : FloatingPointControl;
+
+ auto old = FloatingPointControl.getControlState();
+ FloatingPointControl.setControlState(
+ (old & (-1 - FloatingPointControl.roundingMask)) | FloatingPointControl.roundToZero
+ );
+ x = core.math.rint((x >= 0) ? x + 0.5 : x - 0.5);
+ FloatingPointControl.setControlState(old);
+ return x;
+ }
+ else
+ {
+ return core.stdc.math.roundl(x);
+ }
+}
+
+///
+@safe nothrow @nogc unittest
+{
+ assert(round(4.5) == 5);
+ assert(round(5.4) == 5);
+ assert(round(-4.5) == -5);
+ assert(round(-5.1) == -5);
+}
+
+// assure purity on Posix
+version (Posix)
+{
+ @safe pure nothrow @nogc unittest
+ {
+ assert(round(4.5) == 5);
+ }
+}
+
+/**********************************************
+ * 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 not implemented for Digital Mars C runtime.)
+ */
+long lround(real x) @trusted nothrow @nogc
+{
+ version (CRuntime_DigitalMars)
+ assert(0, "lround not implemented");
+ else
+ return core.stdc.math.llroundl(x);
+}
+
+///
+@safe nothrow @nogc unittest
+{
+ version (CRuntime_DigitalMars) {}
+ else
+ {
+ 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.
+ `pure` on all platforms.
+ */
+real trunc(real x) @trusted nothrow @nogc pure
+{
+ version (InlineAsm_X87_MSVC)
+ {
+ version (X86_64)
+ {
+ 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
+ {
+ 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);
+ }
+}
+
+///
+@safe pure unittest
+{
+ assert(trunc(0.01) == 0);
+ assert(trunc(0.49) == 0);
+ assert(trunc(0.5) == 0);
+ assert(trunc(1.5) == 1);
+}
+
+/*****************************************
+ * 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.
+ */
+pragma(inline, true)
+long rndtol(real x) @nogc @safe pure nothrow { return core.math.rndtol(x); }
+//FIXME
+///ditto
+pragma(inline, true)
+long rndtol(double x) @safe pure nothrow @nogc { return rndtol(cast(real) x); }
+//FIXME
+///ditto
+pragma(inline, true)
+long rndtol(float x) @safe pure nothrow @nogc { return rndtol(cast(real) x); }
+
+///
+@safe unittest
+{
+ assert(rndtol(1.0) == 1L);
+ assert(rndtol(1.2) == 1L);
+ assert(rndtol(1.7) == 2L);
+ assert(rndtol(1.0001) == 1L);
+}
+
+@safe unittest
+{
+ long function(real) prndtol = &rndtol;
+ assert(prndtol != null);
+}
+
+// Helper for floor/ceil
+T floorImpl(T)(const T x) @trusted pure nothrow @nogc
+{
+ import std.math : floatTraits, RealFormat;
+
+ 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;
+
+ // Other kinds of extractors for real formats.
+ static if (F.realFormat == RealFormat.ieeeSingle)
+ int vi;
+ }
+ floatBits y = void;
+ y.rv = x;
+
+ // Find the exponent (power of 2)
+ // Do this by shifting the raw value so that the exponent lies in the low bits,
+ // then mask out the sign bit, and subtract the bias.
+ static if (F.realFormat == RealFormat.ieeeSingle)
+ {
+ int exp = ((y.vi >> (T.mant_dig - 1)) & 0xff) - 0x7f;
+ }
+ 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 ||
+ F.realFormat == RealFormat.ieeeExtended53)
+ {
+ 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;
+ }
+
+ static if (F.realFormat == RealFormat.ieeeSingle)
+ {
+ if (exp < (T.mant_dig - 1))
+ {
+ // Clear all bits representing the fraction part.
+ const uint fraction_mask = F.MANTISSAMASK_INT >> exp;
+
+ if ((y.vi & fraction_mask) != 0)
+ {
+ // If 'x' is negative, then first substract 1.0 from the value.
+ if (y.vi < 0)
+ y.vi += 0x00800000 >> exp;
+ y.vi &= ~fraction_mask;
+ }
+ }
+ }
+ else
+ {
+ static if (F.realFormat == RealFormat.ieeeExtended53)
+ exp = (T.mant_dig + 11 - 1) - exp; // mant_dig is really 64
+ else
+ 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;
+}