Wolfgang Thallergcc16 merge934c3885 months agoHistory
4318 lines · 142.2 KB · D
RawDownload
// Written in the D programming language. /**This is a submodule of $(MREF std, math). It contains several exponential and logarithm functions. Copyright: Copyright The D Language Foundation 2000 - 2011.           D implementations of exp, expm1, exp2, log, log10, log1p, and log2           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 NadlingerSource: $(PHOBOSSRC std/math/exponential.d) Macros:    TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">               <caption>Special Values</caption>               $0</table>    NAN = $(RED NAN)    PLUSMN = &plusmn;    INFIN = &infin;    PLUSMNINF = &plusmn;&infin;    LT = &lt;    GT = &gt; */ module std.math.exponential; import std.traits :  isFloatingPoint, isIntegral, isSigned, isUnsigned, Largest, Unqual; static import core.math;static import core.stdc.math; version (DigitalMars){    version (OSX) { }             // macOS 13 (M1) has issues emulating instruction    else version = INLINE_YL2X;   // x87 has opcodes for these} 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;} version (D_HardFloat){    // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport    version (IeeeFlagsSupport) version = FloatingPointControlSupport;} /** * Compute the value of x $(SUPERSCRIPT n), where n is an integer */Unqual!F pow(F, G)(F x, G n) @nogc @trusted pure nothrowif (isFloatingPoint!(F) && isIntegral!(G)){    import std.traits : Unsigned;     real p = 1.0, v = void;    Unsigned!(Unqual!G) m = n;     if (n < 0)    {        if (n == -1) return 1 / x;         m = cast(typeof(m))(0 - n);        v = p / x;    }    else    {        switch (n)        {        case 0:            return 1.0;        case 1:            return x;        case 2:            return x * x;        default:        }         v = x;    }     while (1)    {        if (m & 1)            p *= v;        m >>= 1;        if (!m)            break;        v *= v;    }    return p;} ///@safe pure nothrow @nogc unittest{    import std.math.operations : feqrel;     assert(pow(2.0, 5) == 32.0);    assert(pow(1.5, 9).feqrel(38.4433) > 16);    assert(pow(real.nan, 2) is real.nan);    assert(pow(real.infinity, 2) == real.infinity);} @safe pure nothrow @nogc unittest{    import std.math.operations : isClose, feqrel;     // Make sure it instantiates and works properly on immutable values and    // with various integer and float types.    immutable real x = 46;    immutable float xf = x;    immutable double xd = x;    immutable uint one = 1;    immutable ushort two = 2;    immutable ubyte three = 3;    immutable ulong eight = 8;     immutable int neg1 = -1;    immutable short neg2 = -2;    immutable byte neg3 = -3;    immutable long neg8 = -8;      assert(pow(x,0) == 1.0);    assert(pow(xd,one) == x);    assert(pow(xf,two) == x * x);    assert(pow(x,three) == x * x * x);    assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x));     assert(pow(x, neg1) == 1 / x);     assert(isClose(pow(xd, neg2), cast(double) (1 / (x * x)), 1e-15));    assert(isClose(pow(xf, neg8), cast(float) (1 / ((x * x) * (x * x) * (x * x) * (x * x))), 1e-15));     assert(feqrel(pow(x, neg3),  1 / (x * x * x)) >= real.mant_dig - 1);} @safe @nogc nothrow unittest{    import std.math.operations : isClose;     assert(isClose(pow(2.0L, 10L), 1024, 1e-18));} // https://issues.dlang.org/show_bug.cgi?id=21601@safe @nogc nothrow pure unittest{    // When reals are large enough the results of pow(b, e) can be    // calculated correctly, if b is of type float or double and e is    // not too large.    static if (real.mant_dig >= 64)    {        // expected result: 3.790e-42        assert(pow(-513645318757045764096.0f, -2) > 0.0);         // expected result: 3.763915357831797e-309        assert(pow(-1.6299717435255677e+154, -2) > 0.0);    }} @safe @nogc nothrow unittest{    import std.math.operations : isClose;    import std.math.traits : isInfinity;     static float f1 = 19100.0f;    static float f2 = 0.000012f;     assert(isClose(pow(f1,9), 3.3829868e+38f));    assert(isInfinity(pow(f1,10)));    assert(pow(f2,9) > 0.0f);    assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));     static double d1 = 21800.0;    static double d2 = 0.000012;     assert(isClose(pow(d1,71), 1.0725339442974e+308));    assert(isInfinity(pow(d1,72)));    assert(pow(d2,65) > 0.0f);    assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));     static if (real.mant_dig == 64) // x87    {        static real r1 = 21950.0L;        static real r2 = 0.000011L;         assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));        assert(isInfinity(pow(r1,1137)));        assert(pow(r2,998) > 0.0L);        assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));    }} @safe @nogc nothrow pure unittest{    import std.math.operations : isClose;     enum f1 = 19100.0f;    enum f2 = 0.000012f;     static assert(isClose(pow(f1,9), 3.3829868e+38f));    static assert(pow(f1,10) > float.max);    static assert(pow(f2,9) > 0.0f);    static assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));     enum d1 = 21800.0;    enum d2 = 0.000012;     static assert(isClose(pow(d1,71), 1.0725339442974e+308));    static assert(pow(d1,72) > double.max);    static assert(pow(d2,65) > 0.0f);    static assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));     static if (real.mant_dig == 64) // x87    {        enum r1 = 21950.0L;        enum r2 = 0.000011L;         static assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));        static assert(pow(r1,1137) > real.max);        static assert(pow(r2,998) > 0.0L);        static assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));    }} /** * Compute the power of two integral numbers. * * Params: *     x = base *     n = exponent * * Returns: *     x raised to the power of n. If n is negative the result is 1 / pow(x, -n), *     which is calculated as integer division with remainder. This may result in *     a division by zero error. * *     If both x and n are 0, the result is 1. * * Throws: *     If x is 0 and n is negative, the result is the same as the result of a *     division by zero. */typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @nogc @safe pure nothrowif (isIntegral!(F) && isIntegral!(G)){    import std.traits : isSigned;     typeof(return) p, v = void;    Unqual!G m = n;     static if (isSigned!(F))    {        if (x == -1) return cast(typeof(return)) (m & 1 ? -1 : 1);    }    static if (isSigned!(G))    {        if (x == 0 && m <= -1) return x / 0;    }    if (x == 1) return 1;    static if (isSigned!(G))    {        if (m < 0) return 0;    }     switch (m)    {    case 0:        p = 1;        break;     case 1:        p = x;        break;     case 2:        p = x * x;        break;     default:        v = x;        p = 1;        while (1)        {            if (m & 1)                p *= v;            m >>= 1;            if (!m)                break;            v *= v;        }        break;    }    return p;} ///@safe pure nothrow @nogc unittest{    assert(pow(2, 3) == 8);    assert(pow(3, 2) == 9);     assert(pow(2, 10) == 1_024);    assert(pow(2, 20) == 1_048_576);    assert(pow(2, 30) == 1_073_741_824);     assert(pow(0, 0) == 1);     assert(pow(1, -5) == 1);    assert(pow(1, -6) == 1);    assert(pow(-1, -5) == -1);    assert(pow(-1, -6) == 1);     assert(pow(-2, 5) == -32);    assert(pow(-2, -5) == 0);    assert(pow(cast(double) -2, -5) == -0.03125);} @safe pure nothrow @nogc unittest{    immutable int one = 1;    immutable byte two = 2;    immutable ubyte three = 3;    immutable short four = 4;    immutable long ten = 10;     assert(pow(two, three) == 8);    assert(pow(two, ten) == 1024);    assert(pow(one, ten) == 1);    assert(pow(ten, four) == 10_000);    assert(pow(four, 10) == 1_048_576);    assert(pow(three, four) == 81);} // https://issues.dlang.org/show_bug.cgi?id=7006@safe pure nothrow @nogc unittest{    assert(pow(5, -1) == 0);    assert(pow(-5, -1) == 0);    assert(pow(5, -2) == 0);    assert(pow(-5, -2) == 0);    assert(pow(-1, int.min) == 1);    assert(pow(-2, int.min) == 0);     assert(pow(4294967290UL,2) == 18446744022169944100UL);    assert(pow(0,uint.max) == 0);} /**Computes integer to floating point powers.*/real pow(I, F)(I x, F y) @nogc @trusted pure nothrowif (isIntegral!I && isFloatingPoint!F){    return pow(cast(real) x, cast(Unqual!F) y);} ///@safe pure nothrow @nogc unittest{    assert(pow(2, 5.0) == 32.0);    assert(pow(7, 3.0) == 343.0);    assert(pow(2, real.nan) is real.nan);    assert(pow(2, real.infinity) == real.infinity);} /** * Calculates x$(SUPERSCRIPT y). * * $(TABLE_SV * $(TR $(TH x) $(TH y) $(TH pow(x, y)) *      $(TH div 0) $(TH invalid?)) * $(TR $(TD anything)      $(TD $(PLUSMN)0.0)                $(TD 1.0) *      $(TD no)        $(TD no) ) * $(TR $(TD |x| $(GT) 1)    $(TD +$(INFIN))                  $(TD +$(INFIN)) *      $(TD no)        $(TD no) ) * $(TR $(TD |x| $(LT) 1)    $(TD +$(INFIN))                  $(TD +0.0) *      $(TD no)        $(TD no) ) * $(TR $(TD |x| $(GT) 1)    $(TD -$(INFIN))                  $(TD +0.0) *      $(TD no)        $(TD no) ) * $(TR $(TD |x| $(LT) 1)    $(TD -$(INFIN))                  $(TD +$(INFIN)) *      $(TD no)        $(TD no) ) * $(TR $(TD +$(INFIN))      $(TD $(GT) 0.0)                  $(TD +$(INFIN)) *      $(TD no)        $(TD no) ) * $(TR $(TD +$(INFIN))      $(TD $(LT) 0.0)                  $(TD +0.0) *      $(TD no)        $(TD no) ) * $(TR $(TD -$(INFIN))      $(TD odd integer $(GT) 0.0)      $(TD -$(INFIN)) *      $(TD no)        $(TD no) ) * $(TR $(TD -$(INFIN))      $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN)) *      $(TD no)        $(TD no)) * $(TR $(TD -$(INFIN))      $(TD odd integer $(LT) 0.0)      $(TD -0.0) *      $(TD no)        $(TD no) ) * $(TR $(TD -$(INFIN))      $(TD $(LT) 0.0, not odd integer) $(TD +0.0) *      $(TD no)        $(TD no) ) * $(TR $(TD $(PLUSMN)1.0)   $(TD $(PLUSMN)$(INFIN))          $(TD -$(NAN)) *      $(TD no)        $(TD yes) ) * $(TR $(TD $(LT) 0.0)      $(TD finite, nonintegral)        $(TD $(NAN)) *      $(TD no)        $(TD yes)) * $(TR $(TD $(PLUSMN)0.0)   $(TD odd integer $(LT) 0.0)      $(TD $(PLUSMNINF)) *      $(TD yes)       $(TD no) ) * $(TR $(TD $(PLUSMN)0.0)   $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN)) *      $(TD yes)       $(TD no)) * $(TR $(TD $(PLUSMN)0.0)   $(TD odd integer $(GT) 0.0)      $(TD $(PLUSMN)0.0) *      $(TD no)        $(TD no) ) * $(TR $(TD $(PLUSMN)0.0)   $(TD $(GT) 0.0, not odd integer) $(TD +0.0) *      $(TD no)        $(TD no) ) * ) */Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrowif (isFloatingPoint!(F) && isFloatingPoint!(G)){    return _powImpl(x, y);} ///@safe pure nothrow @nogc unittest{    import std.math.operations : isClose;     assert(isClose(pow(2.0, 3.0), 8.0));    assert(isClose(pow(1.5, 10.0), 57.6650390625));     // square root of 9    assert(isClose(pow(9.0, 0.5), 3.0));    // 10th root of 1024    assert(isClose(pow(1024.0, 0.1), 2.0));     assert(isClose(pow(-4.0, 3.0), -64.0));     // reciprocal of 4 ^^ 2    assert(isClose(pow(4.0, -2.0), 0.0625));    // reciprocal of (-2) ^^ 3    assert(isClose(pow(-2.0, -3.0), -0.125));     assert(isClose(pow(-2.5, 3.0), -15.625));    // reciprocal of 2.5 ^^ 3    assert(isClose(pow(2.5, -3.0), 0.064));    // reciprocal of (-2.5) ^^ 3    assert(isClose(pow(-2.5, -3.0), -0.064));     // reciprocal of square root of 4    assert(isClose(pow(4.0, -0.5), 0.5));     // per definition    assert(isClose(pow(0.0, 0.0), 1.0));} ///@safe pure nothrow @nogc unittest{    import std.math.operations : isClose;     // the result is a complex number    // which cannot be represented as floating point number    import std.math.traits : isNaN;    assert(isNaN(pow(-2.5, -1.5)));     // use the ^^-operator of std.complex instead    import std.complex : complex;    auto c1 = complex(-2.5, 0.0);    auto c2 = complex(-1.5, 0.0);    auto result = c1 ^^ c2;    // exact result apparently depends on `real` precision => increased tolerance    assert(isClose(result.re, -4.64705438e-17, 2e-4));    assert(isClose(result.im, 2.52982e-1, 2e-4));} @safe pure nothrow @nogc unittest{    import std.math.traits : isNaN;     assert(pow(1.5, real.infinity) == real.infinity);    assert(pow(0.5, real.infinity) == 0.0);    assert(pow(1.5, -real.infinity) == 0.0);    assert(pow(0.5, -real.infinity) == real.infinity);    assert(pow(real.infinity, 1.0) == real.infinity);    assert(pow(real.infinity, -1.0) == 0.0);    assert(pow(real.infinity, real.infinity) == real.infinity);    assert(pow(-real.infinity, 1.0) == -real.infinity);    assert(pow(-real.infinity, 2.0) == real.infinity);    assert(pow(-real.infinity, -1.0) == -0.0);    assert(pow(-real.infinity, -2.0) == 0.0);    assert(isNaN(pow(1.0, real.infinity)));    assert(pow(0.0, -1.0) == real.infinity);    assert(pow(real.nan, 0.0) == 1.0);    assert(isNaN(pow(real.nan, 3.0)));    assert(isNaN(pow(3.0, real.nan)));

Lines 1–500 of 4318.

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