259 lines
7.1 KiB
D
259 lines
7.1 KiB
D
module std.math.traits;
|
||
|
||
import std.traits;
|
||
|
||
enum RealFormat
|
||
{
|
||
ieeeHalf,
|
||
ieeeSingle,
|
||
ieeeDouble,
|
||
ieeeExtended, // x87 80-bit real
|
||
ieeeExtended53, // x87 real rounded to precision of double.
|
||
ibmExtended, // IBM 128-bit extended
|
||
ieeeQuadruple,
|
||
}
|
||
|
||
template floatTraits(T)
|
||
{
|
||
import std.traits : Unqual;
|
||
|
||
// 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 Unqual!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 ulong MANTISSAMASK_LONG = 0x000F_FFFF_FFFF_FFFF;
|
||
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;
|
||
|
||
// For IBM doubledouble the larger magnitude double comes first.
|
||
// It's really a double[2] and arrays don't index differently
|
||
// between little and big-endian targets.
|
||
enum DOUBLEPAIR_MSB = 0;
|
||
enum DOUBLEPAIR_LSB = 1;
|
||
|
||
// The exponent/sign byte is for most significant part.
|
||
version (LittleEndian)
|
||
{
|
||
enum EXPPOS_SHORT = 3;
|
||
enum SIGNPOS_BYTE = 7;
|
||
}
|
||
else
|
||
{
|
||
enum EXPPOS_SHORT = 0;
|
||
enum SIGNPOS_BYTE = 0;
|
||
}
|
||
}
|
||
else
|
||
static assert(false, "No traits support for " ~ T.stringof);
|
||
}
|
||
|
||
bool isNaN(X)(X x) @nogc @trusted pure nothrow
|
||
if (isFloatingPoint!(X))
|
||
{
|
||
version (all)
|
||
{
|
||
return x != x;
|
||
}
|
||
else
|
||
{
|
||
/*
|
||
Code kept for historical context. At least on Intel, the simple test
|
||
x != x uses one dedicated instruction (ucomiss/ucomisd) that runs in one
|
||
cycle. Code for 80- and 128-bits is larger but still smaller than the
|
||
integrals-based solutions below. Future revisions may enable the code
|
||
below conditionally depending on hardware.
|
||
*/
|
||
alias F = floatTraits!(X);
|
||
static if (F.realFormat == RealFormat.ieeeSingle)
|
||
{
|
||
const uint p = *cast(uint *)&x;
|
||
// Sign bit (MSB) is irrelevant so mask it out.
|
||
// Next 8 bits should be all set.
|
||
// At least one bit among the least significant 23 bits should be set.
|
||
return (p & 0x7FFF_FFFF) > 0x7F80_0000;
|
||
}
|
||
else static if (F.realFormat == RealFormat.ieeeDouble)
|
||
{
|
||
const ulong p = *cast(ulong *)&x;
|
||
// Sign bit (MSB) is irrelevant so mask it out.
|
||
// Next 11 bits should be all set.
|
||
// At least one bit among the least significant 52 bits should be set.
|
||
return (p & 0x7FFF_FFFF_FFFF_FFFF) > 0x7FF0_0000_0000_0000;
|
||
}
|
||
else static if (F.realFormat == RealFormat.ieeeExtended ||
|
||
F.realFormat == RealFormat.ieeeExtended53)
|
||
{
|
||
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;
|
||
}
|
||
}
|
||
}
|
||
|
||
bool isInfinity(X)(X x) @nogc @trusted pure nothrow
|
||
if (isFloatingPoint!(X))
|
||
{
|
||
import std.math.traits : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
|
||
|
||
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 ||
|
||
F.realFormat == RealFormat.ieeeExtended53)
|
||
{
|
||
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.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);
|
||
}
|
||
}
|
||
|
||
version (LittleEndian)
|
||
{
|
||
enum MANTISSA_LSB = 0;
|
||
enum MANTISSA_MSB = 1;
|
||
}
|
||
else
|
||
{
|
||
enum MANTISSA_LSB = 1;
|
||
enum MANTISSA_MSB = 0;
|
||
}
|