// Written in the D programming language.
/**
This is a submodule of $(MREF std, math).
It contains several functions for work with floating point numbers.
Copyright: Copyright The D Language Foundation 2000 - 2011.
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/operations.d)
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))
NAN = $(RED NAN)
PLUSMN = ±
INFIN = ∞
LT = <
GT = >
*/
module std.math.operations;
import std.traits : CommonType, isFloatingPoint, isIntegral, Unqual;
// 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
{
import std.math : floatTraits, RealFormat;
alias F = floatTraits!(real);
static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{
// 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
}
if (__ctfe)
{
v = 1; // We use a double in CTFE.
assert(payload >>> 51 == 0,
"Cannot set more than 51 bits of NaN payload in CTFE.");
}
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;
if (__ctfe)
{
v |= 0x7FF0_0000_0000_0000;
return *cast(double*) &v;
}
else 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;
}
}
///
@safe @nogc pure nothrow unittest
{
import std.math.traits : isNaN;
real a = NaN(1_000_000);
assert(isNaN(a));
assert(getNaNPayload(a) == 1_000_000);
}
@system pure nothrow @nogc unittest // not @safe because taking address of local.
{
import std.math : floatTraits, RealFormat;
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
{
import std.math : floatTraits, RealFormat;
// assert(isNaN(x));
alias F = floatTraits!(real);
ulong m = void;
if (__ctfe)
{
double y = x;
m = *cast(ulong*) &y;
// 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.ieeeDouble)
{
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)
{
m = *cast(ulong*)(6+cast(ubyte*)(&x));
}
else
{
m = *cast(ulong*)(2+cast(ubyte*)(&x));
}
m >>= 1; // there's no implicit bit
}
else
{
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;
}
///
@safe @nogc pure nothrow unittest
{
import std.math.traits : isNaN;
real a = NaN(1_000_000);
assert(isNaN(a));
assert(getNaNPayload(a) == 1_000_000);
}
@safe @nogc pure nothrow unittest
{
import std.math.traits : isIdentical, isNaN;
enum real a = NaN(1_000_000);
static assert(isNaN(a));
static assert(getNaNPayload(a) == 1_000_000);
real b = NaN(1_000_000);
assert(isIdentical(b, a));
// The CTFE version of getNaNPayload relies on it being impossible
// for a CTFE-constructed NaN to have more than 51 bits of payload.
enum nanNaN = NaN(getNaNPayload(real.nan));
assert(isIdentical(real.nan, nanNaN));
static if (real.init != real.init)
{
enum initNaN = NaN(getNaNPayload(real.init));
assert(isIdentical(real.init, initNaN));
}
}
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
{
import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
alias F = floatTraits!(real);
static if (F.realFormat != RealFormat.ieeeDouble)
{
if (__ctfe)
{
if (x == -real.infinity)
return -real.max;
if (!(x < real.infinity)) // Infinity or NaN.
return x;
real delta;
// Start with a decent estimate of delta.
if (x <= 0x1.ffffffffffffep+1023 && x >= -double.max)
{
const double d = cast(double) x;
delta = (cast(real) nextUp(d) - cast(real) d) * 0x1p-11L;
while (x + (delta * 0x1p-100L) > x)
delta *= 0x1p-100L;
}
else
{
delta = 0x1p960L;
while (!(x + delta > x) && delta < real.max * 0x1p-100L)
delta *= 0x1p100L;
}
if (x + delta > x)
{
while (x + (delta / 2) > x)
delta /= 2;
}
else
{
do { delta += delta; } while (!(x + delta > x));
}
if (x < 0 && x + delta == 0)
return -0.0L;
return x + delta;
}
}
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 ||
F.realFormat == RealFormat.ieeeExtended53)
{
// For 80-bit reals, the "implied bit" is a nuisance...
ushort *pe = cast(ushort *)&x;
ulong *ps = cast(ulong *)&x;
// EPSILON is 1 for 64-bit, and 2048 for 53-bit precision reals.
enum ulong EPSILON = 2UL ^^ (64 - real.mant_dig);
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 -= EPSILON;
// 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 += EPSILON;
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 s = *cast(ulong *)&x;
if ((s & 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 (s & 0x8000_0000_0000_0000) // Negative number
{
if (s == 0x8000_0000_0000_0000) // it was negative zero
{
s = 0x0000_0000_0000_0001; // change to smallest subnormal
return *cast(double*) &s;
}
--s;
}
else
{ // Positive number
++s;
}
return *cast(double*) &s;
}
/** ditto */
float nextUp(float x) @trusted pure nothrow @nogc
{
uint s = *cast(uint *)&x;
if ((s & 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 (s & 0x8000_0000) // Negative number
{
if (s == 0x8000_0000) // it was negative zero
{
s = 0x0000_0001; // change to smallest subnormal
return *cast(float*) &s;
}
--s;
}
else
{
// Positive number
++s;
}
return *cast(float*) &s;
}
///
@safe @nogc pure nothrow unittest
{
assert(nextUp(1.0 - 1.0e-6).feqrel(0.999999) > 16);
assert(nextUp(1.0 - real.epsilon).feqrel(1.0) > 16);
}
/**
* 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
{
import std.math : floatTraits, RealFormat;
import std.math.traits : isIdentical;
static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
floatTraits!(real).realFormat == RealFormat.ieeeDouble ||
floatTraits!(real).realFormat == RealFormat.ieeeExtended53 ||
floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
{
// Tests for reals
assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC)));
//static 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);
static assert( nextUp(-real.infinity) == -real.max );
static assert( nextUp(-1.0L-real.epsilon) == -1.0 );
static 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) );
static assert( nextUp(-real.min_normal) == -real.min_normal*(1-real.epsilon) );
static assert( nextUp(-real.min_normal*(1-real.epsilon)) == -real.min_normal*(1-2*real.epsilon) );
static assert( -0.0L is nextUp(-real.min_normal*real.epsilon) );
static assert( nextUp(-0.0L) == real.min_normal*real.epsilon );
static assert( nextUp(0.0L) == real.min_normal*real.epsilon );
static assert( nextUp(real.min_normal*(1-real.epsilon)) == real.min_normal );
static 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 );
static assert( nextUp(1.0L) == 1.0 + real.epsilon );
static assert( nextUp(2.0L-real.epsilon) == 2.0 );
static assert( nextUp(real.max) == real.infinity );
static assert( nextUp(real.infinity)==real.infinity );
// ctfe near double.max boundary
static assert(nextUp(nextDown(cast(real) double.max)) == cast(real) double.max);
}
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);
// CTFE
enum double ctfe_n = NaN(0xABC);
//static assert(isIdentical(nextUp(ctfe_n), ctfe_n)); // FIXME: https://issues.dlang.org/show_bug.cgi?id=20197
static assert(nextUp(double.nan) is double.nan);
// negative numbers
static assert( nextUp(-double.infinity) == -double.max );
static assert( nextUp(-1-double.epsilon) == -1.0 );
static assert( nextUp(-2.0) == -2.0 + double.epsilon);
// subnormals and zero
static assert( nextUp(-double.min_normal) == -double.min_normal*(1-double.epsilon) );
static assert( nextUp(-double.min_normal*(1-double.epsilon)) == -double.min_normal*(1-2*double.epsilon) );
static assert( -0.0 is nextUp(-double.min_normal*double.epsilon) );
static assert( nextUp(0.0) == double.min_normal*double.epsilon );
static assert( nextUp(-0.0) == double.min_normal*double.epsilon );
static assert( nextUp(double.min_normal*(1-double.epsilon)) == double.min_normal );
static assert( nextUp(double.min_normal) == double.min_normal*(1+double.epsilon) );
// positive numbers
static assert( nextUp(1.0) == 1.0 + double.epsilon );
static assert( nextUp(2.0-double.epsilon) == 2.0 );
static assert( nextUp(double.max) == double.infinity );
enum float ctfe_fn = NaN(0xABC);
//static assert(isIdentical(nextUp(ctfe_fn), ctfe_fn)); // FIXME: https://issues.dlang.org/show_bug.cgi?id=20197
static assert(nextUp(float.nan) is float.nan);
static assert(nextUp(-float.min_normal) == -float.min_normal*(1-float.epsilon));
static assert(nextUp(1.0f) == 1.0f+float.epsilon);
static assert(nextUp(-0.0f) == float.min_normal*float.epsilon);
static assert(nextUp(float.infinity)==float.infinity);
static assert(nextDown(1.0L+real.epsilon)==1.0);
static assert(nextDown(1.0+double.epsilon)==1.0);
static assert(nextDown(1.0f+float.epsilon)==1.0);
static 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.
* If x or y is a NaN, the result is a NaN.
*
* 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
{
import std.math.traits : isNaN;
if (x == y || isNaN(y))
{
return y;
}
if (isNaN(x))
{
return x;
}
return ((y>x) ? nextUp(x) : nextDown(x));
}
///
@safe pure nothrow @nogc unittest
{
import std.math.traits : isNaN;
float a = 1;
assert(is(typeof(nextafter(a, a)) == float));
assert(nextafter(a, a.infinity) > a);
assert(isNaN(nextafter(a, a.nan)));
assert(isNaN(nextafter(a.nan, a)));
double b = 2;
assert(is(typeof(nextafter(b, b)) == double));
assert(nextafter(b, b.infinity) > b);
assert(isNaN(nextafter(b, b.nan)));
assert(isNaN(nextafter(b.nan, b)));
real c = 3;
assert(is(typeof(nextafter(c, c)) == real));
assert(nextafter(c, c.infinity) > c);
assert(isNaN(nextafter(c, c.nan)));
assert(isNaN(nextafter(c.nan, c)));
}
@safe pure nothrow @nogc unittest
{
import std.math.traits : isNaN, signbit;
// CTFE
enum float a = 1;
static assert(is(typeof(nextafter(a, a)) == float));
static assert(nextafter(a, a.infinity) > a);
static assert(isNaN(nextafter(a, a.nan)));
static assert(isNaN(nextafter(a.nan, a)));
enum double b = 2;
static assert(is(typeof(nextafter(b, b)) == double));
static assert(nextafter(b, b.infinity) > b);
static assert(isNaN(nextafter(b, b.nan)));
static assert(isNaN(nextafter(b.nan, b)));
enum real c = 3;
static assert(is(typeof(nextafter(c, c)) == real));
static assert(nextafter(c, c.infinity) > c);
static assert(isNaN(nextafter(c, c.nan)));
static assert(isNaN(nextafter(c.nan, c)));
enum real negZero = nextafter(+0.0L, -0.0L);
static assert(negZero == -0.0L);
static assert(signbit(negZero));
static assert(nextafter(c, c) == c);
}
//real nexttoward(real x, real y) { return core.stdc.math.nexttowardl(x, y); }
/**
* Returns the positive difference between x and y.
*
* Equivalent to `fmax(x-y, 0)`.
*
* 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) ? +0.0 : x - y;
}
///
@safe pure nothrow @nogc unittest
{
import std.math.traits : isNaN;
assert(fdim(2.0, 0.0) == 2.0);
assert(fdim(-2.0, 0.0) == 0.0);
assert(fdim(real.infinity, 2.0) == real.infinity);
assert(isNaN(fdim(real.nan, 2.0)));
assert(isNaN(fdim(2.0, real.nan)));
assert(isNaN(fdim(real.nan, real.nan)));
}
/**
* Returns the larger of `x` and `y`.
*
* If one of the arguments is a `NaN`, the other is returned.
*
* See_Also: $(REF max, std,algorithm,comparison) is faster because it does not perform the `isNaN` test.
*/
F fmax(F)(const F x, const F y) @safe pure nothrow @nogc
if (__traits(isFloating, F))
{
import std.math.traits : isNaN;
// Do the more predictable test first. Generates 0 branches with ldc and 1 branch with gdc.
// See https://godbolt.org/z/erxrW9
if (isNaN(x)) return y;
return y > x ? y : x;
}
///
@safe pure nothrow @nogc unittest
{
import std.meta : AliasSeq;
static foreach (F; AliasSeq!(float, double, real))
{
assert(fmax(F(0.0), F(2.0)) == 2.0);
assert(fmax(F(-2.0), 0.0) == F(0.0));
assert(fmax(F.infinity, F(2.0)) == F.infinity);
assert(fmax(F.nan, F(2.0)) == F(2.0));
assert(fmax(F(2.0), F.nan) == F(2.0));
}
}
/**
* Returns the smaller of `x` and `y`.
*
* If one of the arguments is a `NaN`, the other is returned.
*
* See_Also: $(REF min, std,algorithm,comparison) is faster because it does not perform the `isNaN` test.
*/
F fmin(F)(const F x, const F y) @safe pure nothrow @nogc
if (__traits(isFloating, F))
{
import std.math.traits : isNaN;
// Do the more predictable test first. Generates 0 branches with ldc and 1 branch with gdc.
// See https://godbolt.org/z/erxrW9
if (isNaN(x)) return y;
return y < x ? y : x;
}
///
@safe pure nothrow @nogc unittest
{
import std.meta : AliasSeq;
static foreach (F; AliasSeq!(float, double, real))
{
assert(fmin(F(0.0), F(2.0)) == 0.0);
assert(fmin(F(-2.0), F(0.0)) == -2.0);
assert(fmin(F.infinity, F(2.0)) == 2.0);
assert(fmin(F.nan, F(2.0)) == 2.0);
assert(fmin(F(2.0), F.nan) == 2.0);
}
}
/**************************************
* Returns (x * y) + z, rounding only once according to the
* current rounding mode.
*
* BUGS: Not currently implemented - rounds twice.
*/
pragma(inline, true)
real fma(real x, real y, real z) @safe pure nothrow @nogc { return (x * y) + z; }
///
@safe pure nothrow @nogc unittest
{
assert(fma(0.0, 2.0, 2.0) == 2.0);
assert(fma(2.0, 2.0, 2.0) == 6.0);
assert(fma(real.infinity, 2.0, 2.0) == real.infinity);
assert(fma(real.nan, 2.0, 2.0) is real.nan);
assert(fma(2.0, 2.0, real.nan) is real.nan);
}
/**************************************
* 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))
{
import std.math : floatTraits, RealFormat;
import core.math : fabs;
/* Public Domain. Author: Don Clugston, 18 Aug 2005.
*/
alias F = floatTraits!(X);
static if (F.realFormat == RealFormat.ieeeSingle
|| F.realFormat == RealFormat.ieeeDouble
|| F.realFormat == RealFormat.ieeeExtended
|| F.realFormat == RealFormat.ieeeExtended53
|| 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;
}
else
{
static assert(false, "Not implemented for this architecture");
}
}
///
@safe pure unittest
{
assert(feqrel(2.0, 2.0) == 53);
assert(feqrel(2.0f, 2.0f) == 24);
assert(feqrel(2.0, double.nan) == 0);
// Test that numbers are within n digits of each
// other by testing if feqrel > n * log2(10)
// five digits
assert(feqrel(2.0, 2.00001) > 16);
// ten digits
assert(feqrel(2.0, 2.00000000001) > 33);
}
@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)();
}
/**
Computes whether a values is approximately equal to a reference value,
admitting a maximum relative difference, and a maximum absolute difference.
Warning:
This template is considered out-dated. It will be removed from
Phobos in 2.106.0. Please use $(LREF isClose) instead. To achieve
a similar behaviour to `approxEqual(a, b)` use
`isClose(a, b, 1e-2, 1e-5)`. In case of comparing to 0.0,
`isClose(a, b, 0.0, eps)` should be used, where `eps`
represents the accepted deviation from 0.0."
Params:
value = Value to compare.
reference = Reference value.
maxRelDiff = Maximum allowable difference relative to `reference`.
Setting to 0.0 disables this check. Defaults to `1e-2`.
maxAbsDiff = Maximum absolute difference. This is mainly usefull
for comparing values to zero. Setting to 0.0 disables this check.
Defaults to `1e-5`.
Returns:
`true` if `value` is approximately equal to `reference` under
either criterium. It is sufficient, when `value ` satisfies
one of the two criteria.
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.
See_Also:
Use $(LREF feqrel) to get the number of equal bits in the mantissa.
*/
deprecated("approxEqual will be removed in 2.106.0. Please use isClose instead.")
bool approxEqual(T, U, V)(T value, U reference, V maxRelDiff = 1e-2, V maxAbsDiff = 1e-5)
{
import core.math : fabs;
import std.range.primitives : empty, front, isInputRange, popFront;
static if (isInputRange!T)
{
static if (isInputRange!U)
{
// Two ranges
for (;; value.popFront(), reference.popFront())
{
if (value.empty) return reference.empty;
if (reference.empty) return value.empty;
if (!approxEqual(value.front, reference.front, maxRelDiff, maxAbsDiff))
return false;
}
}
else static if (isIntegral!U)
{
// convert reference to real
return approxEqual(value, real(reference), maxRelDiff, maxAbsDiff);
}
else
{
// value is range, reference is number
for (; !value.empty; value.popFront())
{
if (!approxEqual(value.front, reference, maxRelDiff, maxAbsDiff))
return false;
}
return true;
}
}
else
{
static if (isInputRange!U)
{
// value is number, reference is range
for (; !reference.empty; reference.popFront())
{
if (!approxEqual(value, reference.front, maxRelDiff, maxAbsDiff))
return false;
}
return true;
}
else static if (isIntegral!T || isIntegral!U)
{
// convert both value and reference to real
return approxEqual(real(value), real(reference), maxRelDiff, maxAbsDiff);
}
else
{
// two numbers
//static assert(is(T : real) && is(U : real));
if (reference == 0)
{
return fabs(value) <= maxAbsDiff;
}
static if (is(typeof(value.infinity)) && is(typeof(reference.infinity)))
{
if (value == value.infinity && reference == reference.infinity ||
value == -value.infinity && reference == -reference.infinity) return true;
}
return fabs((value - reference) / reference) <= maxRelDiff
|| maxAbsDiff != 0 && fabs(value - reference) <= maxAbsDiff;
}
}
}
deprecated @safe pure nothrow unittest
{
assert(approxEqual(1.0, 1.0099));
assert(!approxEqual(1.0, 1.011));
assert(approxEqual(0.00001, 0.0));
assert(!approxEqual(0.00002, 0.0));
assert(approxEqual(3.0, [3, 3.01, 2.99])); // several reference values is strange
assert(approxEqual([3, 3.01, 2.99], 3.0)); // better
float[] arr1 = [ 1.0, 2.0, 3.0 ];
double[] arr2 = [ 1.001, 1.999, 3 ];
assert(approxEqual(arr1, arr2));
}
deprecated @safe pure nothrow unittest
{
// relative comparison depends on reference, make sure proper
// side is used when comparing range to single value. Based on
// https://issues.dlang.org/show_bug.cgi?id=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));
}
deprecated @safe pure nothrow @nogc unittest
{
assert(!approxEqual(0.0,1e-15,1e-9,0.0));
assert(approxEqual(0.0,1e-15,1e-9,1e-9));
assert(!approxEqual(1.0,3.0,0.0,1.0));
assert(approxEqual(1.00000000099,1.0,1e-9,0.0));
assert(!approxEqual(1.0000000011,1.0,1e-9,0.0));
}
deprecated @safe pure nothrow @nogc unittest
{
// maybe unintuitive behavior
assert(approxEqual(1000.0,1010.0));
assert(approxEqual(9_090_000_000.0,9_000_000_000.0));
assert(approxEqual(0.0,1e30,1.0));
assert(approxEqual(0.00001,1e-30));
assert(!approxEqual(-1e-30,1e-30,1e-2,0.0));
}
deprecated @safe pure nothrow @nogc unittest
{
int a = 10;
assert(approxEqual(10, a));
assert(!approxEqual(3, 0));
assert(approxEqual(3, 3));
assert(approxEqual(3.0, 3));
assert(approxEqual(3, 3.0));
assert(approxEqual(0.0,0.0));
assert(approxEqual(-0.0,0.0));
assert(approxEqual(0.0f,0.0));
}
deprecated @safe pure nothrow @nogc unittest
{
real num = real.infinity;
assert(num == real.infinity);
assert(approxEqual(num, real.infinity));
num = -real.infinity;
assert(num == -real.infinity);
assert(approxEqual(num, -real.infinity));
assert(!approxEqual(1,real.nan));
assert(!approxEqual(real.nan,real.max));
assert(!approxEqual(real.nan,real.nan));
}
deprecated @safe pure nothrow unittest
{
assert(!approxEqual([1.0,2.0,3.0],[1.0,2.0]));
assert(!approxEqual([1.0,2.0],[1.0,2.0,3.0]));
assert(approxEqual!(real[],real[])([],[]));
assert(approxEqual(cast(real[])[],cast(real[])[]));
}
/**
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 relative difference.
Setting to 0.0 disables this check. Default depends on the type of
`lhs` and `rhs`: It is approximately half the number of decimal digits of
precision of the smaller type.
maxAbsDiff = Maximum absolute difference. This is mainly usefull
for comparing values to zero. Setting to 0.0 disables this check.
Defaults to `0.0`.
Returns:
`true` if the two items are approximately equal under either criterium.
It is sufficient, when `value ` satisfies one of the two criteria.
If one item is a range, and the other is a single value, then
the result is the logical and-ing of calling `isClose` on
each element of the ranged item against the single item. If
both items are ranges, then `isClose` returns `true` if
and only if the ranges have the same number of elements and if
`isClose` evaluates to `true` for each pair of elements.
See_Also:
Use $(LREF feqrel) to get the number of equal bits in the mantissa.
*/
bool isClose(T, U, V = CommonType!(FloatingPointBaseType!T,FloatingPointBaseType!U))
(T lhs, U rhs, V maxRelDiff = CommonDefaultFor!(T,U), V maxAbsDiff = 0.0)
{
import std.range.primitives : empty, front, isInputRange, popFront;
import std.complex : Complex;
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 (!isClose(lhs.front, rhs.front, maxRelDiff, maxAbsDiff))
return false;
}
}
else
{
// lhs is range, rhs is number
for (; !lhs.empty; lhs.popFront())
{
if (!isClose(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 (!isClose(lhs, rhs.front, maxRelDiff, maxAbsDiff))
return false;
}
return true;
}
else static if (is(T TE == Complex!TE))
{
static if (is(U UE == Complex!UE))
{
// Two complex numbers
return isClose(lhs.re, rhs.re, maxRelDiff, maxAbsDiff)
&& isClose(lhs.im, rhs.im, maxRelDiff, maxAbsDiff);
}
else
{
// lhs is complex, rhs is number
return isClose(lhs.re, rhs, maxRelDiff, maxAbsDiff)
&& isClose(lhs.im, 0.0, maxRelDiff, maxAbsDiff);
}
}
else static if (is(U UE == Complex!UE))
{
// lhs is number, rhs is complex
return isClose(lhs, rhs.re, maxRelDiff, maxAbsDiff)
&& isClose(0.0, rhs.im, maxRelDiff, maxAbsDiff);
}
else
{
// two numbers
if (lhs == rhs) return true;
static if (is(typeof(lhs.infinity)) && is(typeof(rhs.infinity)))
{
if (lhs == lhs.infinity || rhs == rhs.infinity ||
lhs == -lhs.infinity || rhs == -rhs.infinity) return false;
}
import std.math.algebraic : abs;
auto diff = abs(lhs - rhs);
return diff <= maxRelDiff*abs(lhs)
|| diff <= maxRelDiff*abs(rhs)
|| diff <= maxAbsDiff;
}
}
///
@safe pure nothrow @nogc unittest
{
assert(isClose(1.0,0.999_999_999));
assert(isClose(0.001, 0.000_999_999_999));
assert(isClose(1_000_000_000.0,999_999_999.0));
assert(isClose(17.123_456_789, 17.123_456_78));
assert(!isClose(17.123_456_789, 17.123_45));
// use explicit 3rd parameter for less (or more) accuracy
assert(isClose(17.123_456_789, 17.123_45, 1e-6));
assert(!isClose(17.123_456_789, 17.123_45, 1e-7));
// use 4th parameter when comparing close to zero
assert(!isClose(1e-100, 0.0));
assert(isClose(1e-100, 0.0, 0.0, 1e-90));
assert(!isClose(1e-10, -1e-10));
assert(isClose(1e-10, -1e-10, 0.0, 1e-9));
assert(!isClose(1e-300, 1e-298));
assert(isClose(1e-300, 1e-298, 0.0, 1e-200));
// different default limits for different floating point types
assert(isClose(1.0f, 0.999_99f));
assert(!isClose(1.0, 0.999_99));
static if (real.sizeof > double.sizeof)
assert(!isClose(1.0L, 0.999_999_999L));
}
///
@safe pure nothrow unittest
{
assert(isClose([1.0, 2.0, 3.0], [0.999_999_999, 2.000_000_001, 3.0]));
assert(!isClose([1.0, 2.0], [0.999_999_999, 2.000_000_001, 3.0]));
assert(!isClose([1.0, 2.0, 3.0], [0.999_999_999, 2.000_000_001]));
assert(isClose([2.0, 1.999_999_999, 2.000_000_001], 2.0));
assert(isClose(2.0, [2.0, 1.999_999_999, 2.000_000_001]));
}
@safe pure nothrow unittest
{
assert(!isClose([1.0, 2.0, 3.0], [0.999_999_999, 3.0, 3.0]));
assert(!isClose([2.0, 1.999_999, 2.000_000_001], 2.0));
assert(!isClose(2.0, [2.0, 1.999_999_999, 2.000_000_999]));
}
@safe pure nothrow @nogc unittest
{
immutable a = 1.00001f;
const b = 1.000019;
assert(isClose(a,b));
assert(isClose(1.00001f,1.000019f));
assert(isClose(1.00001f,1.000019));
assert(isClose(1.00001,1.000019f));
assert(!isClose(1.00001,1.000019));
real a1 = 1e-300L;
real a2 = a1.nextUp;
assert(isClose(a1,a2));
}
@safe pure nothrow unittest
{
float[] arr1 = [ 1.0, 2.0, 3.0 ];
double[] arr2 = [ 1.00001, 1.99999, 3 ];
assert(isClose(arr1, arr2));
}
@safe pure nothrow @nogc unittest
{
assert(!isClose(1000.0,1010.0));
assert(!isClose(9_090_000_000.0,9_000_000_000.0));
assert(isClose(0.0,1e30,1.0));
assert(!isClose(0.00001,1e-30));
assert(!isClose(-1e-30,1e-30,1e-2,0.0));
}
@safe pure nothrow @nogc unittest
{
assert(!isClose(3, 0));
assert(isClose(3, 3));
assert(isClose(3.0, 3));
assert(isClose(3, 3.0));
assert(isClose(0.0,0.0));
assert(isClose(-0.0,0.0));
assert(isClose(0.0f,0.0));
}
@safe pure nothrow @nogc unittest
{
real num = real.infinity;
assert(num == real.infinity);
assert(isClose(num, real.infinity));
num = -real.infinity;
assert(num == -real.infinity);
assert(isClose(num, -real.infinity));
assert(!isClose(1,real.nan));
assert(!isClose(real.nan,real.max));
assert(!isClose(real.nan,real.nan));
}
@safe pure nothrow @nogc unittest
{
assert(isClose!(real[],real[],real)([],[]));
assert(isClose(cast(real[])[],cast(real[])[]));
}
@safe pure nothrow @nogc unittest
{
import std.conv : to;
float f = 31.79f;
double d = 31.79;
double f2d = f.to!double;
assert(isClose(f,f2d));
assert(!isClose(d,f2d));
}
@safe pure nothrow @nogc unittest
{
import std.conv : to;
double d = 31.79;
float f = d.to!float;
double f2d = f.to!double;
assert(isClose(f,f2d));
assert(!isClose(d,f2d));
assert(isClose(d,f2d,1e-4));
}
package(std.math) template CommonDefaultFor(T,U)
{
import std.algorithm.comparison : min;
alias baseT = FloatingPointBaseType!T;
alias baseU = FloatingPointBaseType!U;
enum CommonType!(baseT, baseU) CommonDefaultFor = 10.0L ^^ -((min(baseT.dig, baseU.dig) + 1) / 2 + 1);
}
private template FloatingPointBaseType(T)
{
import std.range.primitives : ElementType;
static if (isFloatingPoint!T)
{
alias FloatingPointBaseType = Unqual!T;
}
else static if (isFloatingPoint!(ElementType!(Unqual!T)))
{
alias FloatingPointBaseType = Unqual!(ElementType!(Unqual!T));
}
else
{
alias FloatingPointBaseType = real;
}
}
/***********************************
* 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 `x` precedes `y` in the order specified above;
* 0 if `x` and `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)
{
import std.math : floatTraits, RealFormat;
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
import std.math.traits : signbit, isNaN;
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;
static 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);
}
}}
}
package(std): // not yet public
struct FloatingPointBitpattern(T)
if (isFloatingPoint!T)
{
static if (T.mant_dig <= 64)
{
ulong mantissa;
}
else
{
ulong mantissa_lsb;
ulong mantissa_msb;
}
int exponent;
bool negative;
}
FloatingPointBitpattern!T extractBitpattern(T)(const(T) value) @trusted
if (isFloatingPoint!T)
{
import std.math : floatTraits, RealFormat;
T val = value;
FloatingPointBitpattern!T ret;
alias F = floatTraits!T;
static if (F.realFormat == RealFormat.ieeeExtended)
{
if (__ctfe)
{
import core.math : fabs, ldexp;
import std.math.rounding : floor;
import std.math.traits : isInfinity, isNaN, signbit;
import std.math.exponential : log2;
if (isNaN(val) || isInfinity(val))
ret.exponent = 32767;
else if (fabs(val) < real.min_normal)
ret.exponent = 0;
else if (fabs(val) >= nextUp(real.max / 2))
ret.exponent = 32766;
else
ret.exponent = cast(int) (val.fabs.log2.floor() + 16383);
if (ret.exponent == 32767)
{
// NaN or infinity
ret.mantissa = isNaN(val) ? ((1L << 63) - 1) : 0;
}
else
{
auto delta = 16382 + 64 // bias + bits of ulong
- (ret.exponent == 0 ? 1 : ret.exponent); // -1 in case of subnormals
val = ldexp(val, delta); // val *= 2^^delta
ulong tmp = cast(ulong) fabs(val);
if (ret.exponent != 32767 && ret.exponent > 0 && tmp <= ulong.max / 2)
{
// correction, due to log2(val) being rounded up:
ret.exponent--;
val *= 2;
tmp = cast(ulong) fabs(val);
}
ret.mantissa = tmp & long.max;
}
ret.negative = (signbit(val) == 1);
}
else
{
ushort* vs = cast(ushort*) &val;
ret.mantissa = (cast(ulong*) vs)[0] & long.max;
ret.exponent = vs[4] & short.max;
ret.negative = (vs[4] >> 15) & 1;
}
}
else
{
static if (F.realFormat == RealFormat.ieeeSingle)
{
ulong ival = *cast(uint*) &val;
}
else static if (F.realFormat == RealFormat.ieeeDouble)
{
ulong ival = *cast(ulong*) &val;
}
else
{
static assert(false, "Floating point type `" ~ F.realFormat ~ "` not supported.");
}
import std.math.exponential : log2;
enum log2_max_exp = cast(int) log2(T.max_exp);
ret.mantissa = ival & ((1L << (T.mant_dig - 1)) - 1);
ret.exponent = (ival >> (T.mant_dig - 1)) & ((1L << (log2_max_exp + 1)) - 1);
ret.negative = (ival >> (T.mant_dig + log2_max_exp)) & 1;
}
// add leading 1 for normalized values and correct exponent for denormalied values
if (ret.exponent != 0 && ret.exponent != 2 * T.max_exp - 1)
ret.mantissa |= 1L << (T.mant_dig - 1);
else if (ret.exponent == 0)
ret.exponent = 1;
ret.exponent -= T.max_exp - 1;
return ret;
}
@safe pure unittest
{
float f = 1.0f;
auto bp = extractBitpattern(f);
assert(bp.mantissa == 0x80_0000);
assert(bp.exponent == 0);
assert(bp.negative == false);
f = float.max;
bp = extractBitpattern(f);
assert(bp.mantissa == 0xff_ffff);
assert(bp.exponent == 127);
assert(bp.negative == false);
f = -1.5432e-17f;
bp = extractBitpattern(f);
assert(bp.mantissa == 0x8e_55c8);
assert(bp.exponent == -56);
assert(bp.negative == true);
// using double literal due to https://issues.dlang.org/show_bug.cgi?id=20361
f = 2.3822073893521890206e-44;
bp = extractBitpattern(f);
assert(bp.mantissa == 0x00_0011);
assert(bp.exponent == -126);
assert(bp.negative == false);
f = -float.infinity;
bp = extractBitpattern(f);
assert(bp.mantissa == 0);
assert(bp.exponent == 128);
assert(bp.negative == true);
f = float.nan;
bp = extractBitpattern(f);
assert(bp.mantissa != 0); // we don't guarantee payloads
assert(bp.exponent == 128);
assert(bp.negative == false);
}
@safe pure unittest
{
double d = 1.0;
auto bp = extractBitpattern(d);
assert(bp.mantissa == 0x10_0000_0000_0000L);
assert(bp.exponent == 0);
assert(bp.negative == false);
d = double.max;
bp = extractBitpattern(d);
assert(bp.mantissa == 0x1f_ffff_ffff_ffffL);
assert(bp.exponent == 1023);
assert(bp.negative == false);
d = -1.5432e-222;
bp = extractBitpattern(d);
assert(bp.mantissa == 0x11_d9b6_a401_3b04L);
assert(bp.exponent == -737);
assert(bp.negative == true);
d = 0.0.nextUp;
bp = extractBitpattern(d);
assert(bp.mantissa == 0x00_0000_0000_0001L);
assert(bp.exponent == -1022);
assert(bp.negative == false);
d = -double.infinity;
bp = extractBitpattern(d);
assert(bp.mantissa == 0);
assert(bp.exponent == 1024);
assert(bp.negative == true);
d = double.nan;
bp = extractBitpattern(d);
assert(bp.mantissa != 0); // we don't guarantee payloads
assert(bp.exponent == 1024);
assert(bp.negative == false);
}
@safe pure unittest
{
import std.math : floatTraits, RealFormat;
alias F = floatTraits!real;
static if (F.realFormat == RealFormat.ieeeExtended)
{
real r = 1.0L;
auto bp = extractBitpattern(r);
assert(bp.mantissa == 0x8000_0000_0000_0000L);
assert(bp.exponent == 0);
assert(bp.negative == false);
r = real.max;
bp = extractBitpattern(r);
assert(bp.mantissa == 0xffff_ffff_ffff_ffffL);
assert(bp.exponent == 16383);
assert(bp.negative == false);
r = -1.5432e-3333L;
bp = extractBitpattern(r);
assert(bp.mantissa == 0xc768_a2c7_a616_cc22L);
assert(bp.exponent == -11072);
assert(bp.negative == true);
r = 0.0L.nextUp;
bp = extractBitpattern(r);
assert(bp.mantissa == 0x0000_0000_0000_0001L);
assert(bp.exponent == -16382);
assert(bp.negative == false);
r = -float.infinity;
bp = extractBitpattern(r);
assert(bp.mantissa == 0);
assert(bp.exponent == 16384);
assert(bp.negative == true);
r = float.nan;
bp = extractBitpattern(r);
assert(bp.mantissa != 0); // we don't guarantee payloads
assert(bp.exponent == 16384);
assert(bp.negative == false);
r = nextDown(0x1p+16383L);
bp = extractBitpattern(r);
assert(bp.mantissa == 0xffff_ffff_ffff_ffffL);
assert(bp.exponent == 16382);
assert(bp.negative == false);
}
}
@safe pure unittest
{
import std.math : floatTraits, RealFormat;
import std.math.exponential : log2;
alias F = floatTraits!real;
// log2 is broken for x87-reals on some computers in CTFE
// the following test excludes these computers from the test
// (issue 21757)
enum test = cast(int) log2(3.05e2312L);
static if (F.realFormat == RealFormat.ieeeExtended && test == 7681)
{
enum r1 = 1.0L;
enum bp1 = extractBitpattern(r1);
static assert(bp1.mantissa == 0x8000_0000_0000_0000L);
static assert(bp1.exponent == 0);
static assert(bp1.negative == false);
enum r2 = real.max;
enum bp2 = extractBitpattern(r2);
static assert(bp2.mantissa == 0xffff_ffff_ffff_ffffL);
static assert(bp2.exponent == 16383);
static assert(bp2.negative == false);
enum r3 = -1.5432e-3333L;
enum bp3 = extractBitpattern(r3);
static assert(bp3.mantissa == 0xc768_a2c7_a616_cc22L);
static assert(bp3.exponent == -11072);
static assert(bp3.negative == true);
enum r4 = 0.0L.nextUp;
enum bp4 = extractBitpattern(r4);
static assert(bp4.mantissa == 0x0000_0000_0000_0001L);
static assert(bp4.exponent == -16382);
static assert(bp4.negative == false);
enum r5 = -real.infinity;
enum bp5 = extractBitpattern(r5);
static assert(bp5.mantissa == 0);
static assert(bp5.exponent == 16384);
static assert(bp5.negative == true);
enum r6 = real.nan;
enum bp6 = extractBitpattern(r6);
static assert(bp6.mantissa != 0); // we don't guarantee payloads
static assert(bp6.exponent == 16384);
static assert(bp6.negative == false);
enum r7 = nextDown(0x1p+16383L);
enum bp7 = extractBitpattern(r7);
static assert(bp7.mantissa == 0xffff_ffff_ffff_ffffL);
static assert(bp7.exponent == 16382);
static assert(bp7.negative == false);
}
}