// 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 NadlingerSource: $(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 = &plusmn;    INFIN = &infin;    LT = &lt;    GT = &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.traits : 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.traits : 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 @nogcin{    // Precondition: Input must be NaN    import std.math.traits : isNaN;    assert(isNaN(x), "getNaNPayload called on a non-NaN value");}do{    import std.math.traits : 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.traits : 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  )

Lines 1–500 of 2011. Load more