Wolfgang Thallergcc16 merge934c3885 months agoHistory
4201 lines · 129.6 KB · D
RawDownload
// Written in the D programming language. /**This module is a port of a growing fragment of the $(D_PARAM numeric)header in Alexander Stepanov's $(LINK2 https://en.wikipedia.org/wiki/Standard_Template_Library,Standard Template Library), with a few additions. Macros:Copyright: Copyright Andrei Alexandrescu 2008 - 2009.License:   $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).Authors:   $(HTTP erdani.org, Andrei Alexandrescu),                   Don Clugston, Robert Jacques, Ilya YaroshenkoSource:    $(PHOBOSSRC std/numeric.d)*//*         Copyright Andrei Alexandrescu 2008 - 2009.Distributed under the Boost Software License, Version 1.0.   (See accompanying file LICENSE_1_0.txt or copy at         http://www.boost.org/LICENSE_1_0.txt)*/module std.numeric; import std.complex;import std.math;import core.math : fabs, ldexp, sin, sqrt;import std.range.primitives;import std.traits;import std.typecons; /// Format flags for CustomFloat.public enum CustomFloatFlags{    /// Adds a sign bit to allow for signed numbers.    signed = 1,     /**     * Store values in normalized form by default. The actual precision of the     * significand is extended by 1 bit by assuming an implicit leading bit of 1     * instead of 0. i.e. `1.nnnn` instead of `0.nnnn`.     * True for all $(LINK2 https://en.wikipedia.org/wiki/IEEE_floating_point, IEEE754) types     */    storeNormalized = 2,     /**     * Stores the significand in $(LINK2 https://en.wikipedia.org/wiki/IEEE_754-1985#Denormalized_numbers,     * IEEE754 denormalized) form when the exponent is 0. Required to express the value 0.     */    allowDenorm = 4,     /**      * Allows the storage of $(LINK2 https://en.wikipedia.org/wiki/IEEE_754-1985#Positive_and_negative_infinity,      * IEEE754 _infinity) values.      */    infinity = 8,     /// Allows the storage of $(LINK2 https://en.wikipedia.org/wiki/NaN, IEEE754 Not a Number) values.    nan = 16,     /**     * If set, select an exponent bias such that max_exp = 1.     * i.e. so that the maximum value is >= 1.0 and < 2.0.     * Ignored if the exponent bias is manually specified.     */    probability = 32,     /// If set, unsigned custom floats are assumed to be negative.    negativeUnsigned = 64,     /**If set, 0 is the only allowed $(LINK2 https://en.wikipedia.org/wiki/IEEE_754-1985#Denormalized_numbers,     * IEEE754 denormalized) number.     * Requires allowDenorm and storeNormalized.     */    allowDenormZeroOnly = 128 | allowDenorm | storeNormalized,     /// Include _all of the $(LINK2 https://en.wikipedia.org/wiki/IEEE_floating_point, IEEE754) options.    ieee = signed | storeNormalized | allowDenorm | infinity | nan ,     /// Include none of the above options.    none = 0} private enum isIEEEQuadruple = floatTraits!real.realFormat == RealFormat.ieeeQuadruple; private template CustomFloatParams(uint bits){    enum CustomFloatFlags flags = CustomFloatFlags.ieee                ^ ((bits == 80 && !isIEEEQuadruple) ? CustomFloatFlags.storeNormalized : CustomFloatFlags.none);    static if (bits ==  8) alias CustomFloatParams = CustomFloatParams!( 4,  3, flags);    static if (bits == 16) alias CustomFloatParams = CustomFloatParams!(10,  5, flags);    static if (bits == 32) alias CustomFloatParams = CustomFloatParams!(23,  8, flags);    static if (bits == 64) alias CustomFloatParams = CustomFloatParams!(52, 11, flags);    static if (bits == 80) alias CustomFloatParams = CustomFloatParams!(64, 15, flags);} private template CustomFloatParams(uint precision, uint exponentWidth, CustomFloatFlags flags){    import std.meta : AliasSeq;    alias CustomFloatParams =        AliasSeq!(            precision,            exponentWidth,            flags,            (1 << (exponentWidth - ((flags & flags.probability) == 0)))             - ((flags & (flags.nan | flags.infinity)) != 0) - ((flags & flags.probability) != 0)        ); // ((flags & CustomFloatFlags.probability) == 0)} /** * Allows user code to define custom floating-point formats. These formats are * for storage only; all operations on them are performed by first implicitly * extracting them to `real` first. After the operation is completed the * result can be stored in a custom floating-point value via assignment. */template CustomFloat(uint bits)if (bits == 8 || bits == 16 || bits == 32 || bits == 64 || bits == 80){    alias CustomFloat = CustomFloat!(CustomFloatParams!(bits));} /// dittotemplate CustomFloat(uint precision, uint exponentWidth, CustomFloatFlags flags = CustomFloatFlags.ieee)if (((flags & flags.signed) + precision + exponentWidth) % 8 == 0 && precision + exponentWidth > 0){    alias CustomFloat = CustomFloat!(CustomFloatParams!(precision, exponentWidth, flags));} ///@safe unittest{    import std.math.trigonometry : sin, cos;     // Define a 16-bit floating point values    CustomFloat!16                                x;     // Using the number of bits    CustomFloat!(10, 5)                           y;     // Using the precision and exponent width    CustomFloat!(10, 5,CustomFloatFlags.ieee)     z;     // Using the precision, exponent width and format flags    CustomFloat!(10, 5,CustomFloatFlags.ieee, 15) w;     // Using the precision, exponent width, format flags and exponent offset bias     // Use the 16-bit floats mostly like normal numbers    w = x*y - 1;     // Functions calls require conversion    z = sin(+x)           + cos(+y);                     // Use unary plus to concisely convert to a real    z = sin(x.get!float)  + cos(y.get!float);            // Or use get!T    z = sin(cast(float) x) + cos(cast(float) y);           // Or use cast(T) to explicitly convert     // Define a 8-bit custom float for storing probabilities    alias Probability = CustomFloat!(4, 4, CustomFloatFlags.ieee^CustomFloatFlags.probability^CustomFloatFlags.signed );    auto p = Probability(0.5);} // Facilitate converting numeric types to custom floatprivate union ToBinary(F)if (is(typeof(CustomFloatParams!(F.sizeof*8))) || is(F == real)){    F set;     // If on Linux or Mac, where 80-bit reals are padded, ignore the    // padding.    import std.algorithm.comparison : min;    CustomFloat!(CustomFloatParams!(min(F.sizeof*8, 80))) get;     // Convert F to the correct binary type.    static typeof(get) opCall(F value)    {        ToBinary r;        r.set = value;        return r.get;    }    alias get this;} /// dittostruct CustomFloat(uint             precision,  // fraction bits (23 for float)                   uint             exponentWidth,  // exponent bits (8 for float)  Exponent width                   CustomFloatFlags flags,                   uint             bias)if (isCorrectCustomFloat(precision, exponentWidth, flags)){    import std.bitmanip : bitfields;    import std.meta : staticIndexOf;private:    // get the correct unsigned bitfield type to support > 32 bits    template uType(uint bits)    {        static if (bits <= size_t.sizeof*8)  alias uType = size_t;        else                                alias uType = ulong ;    }     // get the correct signed   bitfield type to support > 32 bits    template sType(uint bits)    {        static if (bits <= ptrdiff_t.sizeof*8-1) alias sType = ptrdiff_t;        else                                    alias sType = long;    }     alias T_sig = uType!precision;    alias T_exp = uType!exponentWidth;    alias T_signed_exp = sType!exponentWidth;     alias Flags = CustomFloatFlags;     // Perform IEEE rounding with round to nearest detection    void roundedShift(T,U)(ref T sig, U shift)    {        if (shift >= T.sizeof*8)        {            // avoid illegal shift            sig = 0;        }        else if (sig << (T.sizeof*8 - shift) == cast(T) 1uL << (T.sizeof*8 - 1))        {            // round to even            sig >>= shift;            sig  += sig & 1;        }        else        {            sig >>= shift - 1;            sig  += sig & 1;            // Perform standard rounding            sig >>= 1;        }    }     // Convert the current value to signed exponent, normalized form    void toNormalized(T,U)(ref T sig, ref U exp) const    {        sig = significand;        auto shift = (T.sizeof*8) - precision;        exp = exponent;        static if (flags&(Flags.infinity|Flags.nan))        {            // Handle inf or nan            if (exp == exponent_max)            {                exp = exp.max;                sig <<= shift;                static if (flags&Flags.storeNormalized)                {                    // Save inf/nan in denormalized format                    sig >>= 1;                    sig  += cast(T) 1uL << (T.sizeof*8 - 1);                }                return;            }        }        if ((~flags&Flags.storeNormalized) ||            // Convert denormalized form to normalized form            ((flags&Flags.allowDenorm) && exp == 0))        {            if (sig > 0)            {                import core.bitop : bsr;                auto shift2 = precision - bsr(sig);                exp  -= shift2-1;                shift += shift2;            }            else                                // value = 0.0            {                exp = exp.min;                return;            }        }        sig <<= shift;        exp -= bias;    }     // Set the current value from signed exponent, normalized form    void fromNormalized(T,U)(ref T sig, ref U exp)    {        auto shift = (T.sizeof*8) - precision;        if (exp == exp.max)        {            // infinity or nan            exp = exponent_max;            static if (flags & Flags.storeNormalized)                sig <<= 1;             // convert back to normalized form            static if (~flags & Flags.infinity)                // No infinity support?                assert(sig != 0, "Infinity floating point value assigned to a "                        ~ typeof(this).stringof ~ " (no infinity support).");             static if (~flags & Flags.nan)  // No NaN support?                assert(sig == 0, "NaN floating point value assigned to a " ~                        typeof(this).stringof ~ " (no nan support).");            sig >>= shift;            return;        }        if (exp == exp.min)     // 0.0        {             exp = 0;             sig = 0;             return;        }         exp += bias;        if (exp <= 0)        {            static if ((flags&Flags.allowDenorm) ||                       // Convert from normalized form to denormalized                       (~flags&Flags.storeNormalized))            {                shift += -exp;                roundedShift(sig,1);                sig   += cast(T) 1uL << (T.sizeof*8 - 1);                // Add the leading 1                exp    = 0;            }            else                assert((flags&Flags.storeNormalized) && exp == 0,                    "Underflow occured assigning to a " ~                    typeof(this).stringof ~ " (no denormal support).");        }        else        {            static if (~flags&Flags.storeNormalized)            {                // Convert from normalized form to denormalized                roundedShift(sig,1);                sig  += cast(T) 1uL << (T.sizeof*8 - 1);                // Add the leading 1            }        }         if (shift > 0)            roundedShift(sig,shift);        if (sig > significand_max)        {            // handle significand overflow (should only be 1 bit)            static if (~flags&Flags.storeNormalized)            {                sig >>= 1;            }            else                sig &= significand_max;            exp++;        }        static if ((flags&Flags.allowDenormZeroOnly)==Flags.allowDenormZeroOnly)        {            // disallow non-zero denormals            if (exp == 0)            {                sig <<= 1;                if (sig > significand_max && (sig&significand_max) > 0)                    // Check and round to even                    exp++;                sig = 0;            }        }         if (exp >= exponent_max)        {            static if (flags&(Flags.infinity|Flags.nan))            {                sig         = 0;                exp         = exponent_max;                static if (~flags&(Flags.infinity))                    assert(0, "Overflow occured assigning to a " ~                        typeof(this).stringof ~ " (no infinity support).");            }            else                assert(exp == exponent_max, "Overflow occured assigning to a "                    ~ typeof(this).stringof ~ " (no infinity support).");        }    } public:    static if (precision == 64) // CustomFloat!80 support hack    {        static if (isIEEEQuadruple)        {        // Only use highest 64 significand bits from 112 explicitly stored        align (1):            enum ulong significand_max = ulong.max;            version (LittleEndian)            {                private ubyte[6] _padding; // 48-bit of padding                ulong significand;                mixin(bitfields!(                    T_exp , "exponent", exponentWidth,                    bool  , "sign"    , flags & flags.signed ));            }            else            {                mixin(bitfields!(                    T_exp , "exponent", exponentWidth,                    bool  , "sign"    , flags & flags.signed ));                ulong significand;                private ubyte[6] _padding; // 48-bit of padding            }        }        else        {            ulong significand;            enum ulong significand_max = ulong.max;            mixin(bitfields!(                T_exp , "exponent", exponentWidth,                bool  , "sign"    , flags & flags.signed ));        }    }    else    {        mixin(bitfields!(            T_sig, "significand", precision,            T_exp, "exponent"   , exponentWidth,            bool , "sign"       , flags & flags.signed ));    }     /// Returns: infinity value    static if (flags & Flags.infinity)        static @property CustomFloat infinity()        {            CustomFloat value;            static if (flags & Flags.signed)                value.sign          = 0;            value.significand   = 0;            value.exponent      = exponent_max;            return value;        }     /// Returns: NaN value    static if (flags & Flags.nan)        static @property CustomFloat nan()        {            CustomFloat value;            static if (flags & Flags.signed)                value.sign          = 0;            value.significand   = cast(typeof(significand_max)) 1L << (precision-1);            value.exponent      = exponent_max;            return value;        }     /// Returns: number of decimal digits of precision    static @property size_t dig()    {        auto shiftcnt = precision - ((flags&Flags.storeNormalized) == 0);        return shiftcnt == 64 ? 19 : cast(size_t) log10(real(1uL << shiftcnt));    }     /// Returns: smallest increment to the value 1    static @property CustomFloat epsilon()    {        CustomFloat one = CustomFloat(1);        CustomFloat onePlusEpsilon = one;        onePlusEpsilon.significand = onePlusEpsilon.significand | 1; // |= does not work here         return CustomFloat(onePlusEpsilon - one);    }     /// the number of bits in mantissa    enum mant_dig = precision + ((flags&Flags.storeNormalized) != 0);     /// Returns: maximum int value such that 10<sup>max_10_exp</sup> is representable    static @property int max_10_exp(){ return cast(int) log10( +max ); }     /// maximum int value such that 2<sup>max_exp-1</sup> is representable    enum max_exp = exponent_max - bias - ((flags & (Flags.infinity | Flags.nan)) != 0) + 1;     /// Returns: minimum int value such that 10<sup>min_10_exp</sup> is representable    static @property int min_10_exp(){ return cast(int) log10( +min_normal ); }     /// minimum int value such that 2<sup>min_exp-1</sup> is representable as a normalized value    enum min_exp = cast(T_signed_exp) -(cast(long) bias) + 1 + ((flags & Flags.allowDenorm) != 0);     /// Returns: largest representable value that's not infinity    static @property CustomFloat max()    {        CustomFloat value;        static if (flags & Flags.signed)            value.sign        = 0;        value.exponent    = exponent_max - ((flags&(flags.infinity|flags.nan)) != 0);        value.significand = significand_max;        return value;    }     /// Returns: smallest representable normalized value that's not 0    static @property CustomFloat min_normal()    {        CustomFloat value;        static if (flags & Flags.signed)            value.sign = 0;        value.exponent = (flags & Flags.allowDenorm) != 0;        static if (flags & Flags.storeNormalized)            value.significand = 0;        else            value.significand = cast(T_sig) 1uL << (precision - 1);        return value;    }     /// Returns: real part    @property CustomFloat re() const { return this; }     /// Returns: imaginary part    static @property CustomFloat im() { return CustomFloat(0.0f); }     /// Initialize from any `real` compatible type.    this(F)(F input)    if (__traits(compiles, cast(real) input ))

Lines 1–500 of 4201.

Keyboard shortcuts

General

Show keyboard shortcuts
?
Focus the current search
/

Issues and pull requests

Open filters
f

Issues, pull requests and Inbox

Select next item
j↓
Select previous item
k↑
Open selected item
Enter
Dismiss filters or clear selection
Esc

Repository code

Go to file
t