Wolfgang Thallergcc 15.1.0ada185d1 year agoHistory
1424 lines · 40.2 KB · D
RawDownload
// Written in the D programming language. /**This is a submodule of $(MREF std, math). It contains several trigonometric functions. Copyright: Copyright The D Language Foundation 2000 - 2011.           D implementations of tan, atan, and atan2 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/trigonometry.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))    TH3 = $(TR $(TH $1) $(TH $2) $(TH $3))    TD3 = $(TR $(TD $1) $(TD $2) $(TD $3))    TABLE_DOMRG = <table border="1" cellpadding="4" cellspacing="0">                  $(SVH Domain X, Range Y)                  $(SV $1, $2)                  </table>    DOMAIN=$1    RANGE=$1    POWER = $1<sup>$2</sup>    NAN = $(RED NAN)    PLUSMN = &plusmn;    INFIN = &infin;    PLUSMNINF = &plusmn;&infin; */ module std.math.trigonometry; static import core.math; 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;} /*********************************** * Returns cosine of x. x is in radians. * *      $(TABLE_SV *      $(TR $(TH x)                 $(TH cos(x)) $(TH invalid?)) *      $(TR $(TD $(NAN))            $(TD $(NAN)) $(TD yes)     ) *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)     ) *      ) * Bugs: *      Results are undefined if |x| >= $(POWER 2,64). */pragma(inline, true)real cos(real x) @safe pure nothrow @nogc { return core.math.cos(x); }///dittopragma(inline, true)double cos(double x) @safe pure nothrow @nogc { return core.math.cos(x); }///dittopragma(inline, true)float cos(float x) @safe pure nothrow @nogc { return core.math.cos(x); } ///@safe unittest{    import std.math.operations : isClose;     assert(cos(0.0) == 1.0);    assert(cos(1.0).isClose(0.5403023059));    assert(cos(3.0).isClose(-0.9899924966));} @safe unittest{    real function(real) pcos = &cos;    assert(pcos != null);} @safe pure nothrow @nogc unittest{    import std.math.algebraic : fabs;     float f = cos(-2.0f);    assert(fabs(f - -0.416147f) < .00001);     double d = cos(-2.0);    assert(fabs(d - -0.416147f) < .00001);     real r = cos(-2.0L);    assert(fabs(r - -0.416147f) < .00001);} /*********************************** * Returns $(HTTP en.wikipedia.org/wiki/Sine, sine) of x. x is in $(HTTP en.wikipedia.org/wiki/Radian, radians). * *      $(TABLE_SV *      $(TH3 x           ,  sin(x)      ,  invalid?) *      $(TD3 $(NAN)      ,  $(NAN)      ,  yes     ) *      $(TD3 $(PLUSMN)0.0,  $(PLUSMN)0.0,  no      ) *      $(TD3 $(PLUSMNINF),  $(NAN)      ,  yes     ) *      ) * * Params: *      x = angle in radians (not degrees) * Returns: *      sine of x * See_Also: *      $(MYREF cos), $(MYREF tan), $(MYREF asin) * Bugs: *      Results are undefined if |x| >= $(POWER 2,64). */pragma(inline, true)real sin(real x) @safe pure nothrow @nogc { return core.math.sin(x); }///dittopragma(inline, true)double sin(double x) @safe pure nothrow @nogc { return core.math.sin(x); }///dittopragma(inline, true)float sin(float x) @safe pure nothrow @nogc { return core.math.sin(x); } ///@safe unittest{    import std.math.constants : PI;    import std.stdio : writefln;     void someFunc()    {      real x = 30.0;      auto result = sin(x * (PI / 180)); // convert degrees to radians      writefln("The sine of %s degrees is %s", x, result);    }} @safe unittest{    real function(real) psin = &sin;    assert(psin != null);} @safe pure nothrow @nogc unittest{    import std.math.algebraic : fabs;     float f = sin(-2.0f);    assert(fabs(f - -0.909297f) < .00001);     double d = sin(-2.0);    assert(fabs(d - -0.909297f) < .00001);     real r = sin(-2.0L);    assert(fabs(r - -0.909297f) < .00001);} /**************************************************************************** * Returns tangent of x. x is in radians. * *      $(TABLE_SV *      $(TR $(TH x)             $(TH tan(x))       $(TH invalid?)) *      $(TR $(TD $(NAN))        $(TD $(NAN))       $(TD yes)) *      $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) $(TD no)) *      $(TR $(TD $(PLUSMNINF))  $(TD $(NAN))       $(TD yes)) *      ) */pragma(inline, true)real tan(real x) @safe pure nothrow @nogc{    version (InlineAsm_X87)    {        if (!__ctfe)            return tanAsm(x);    }    return tanImpl(x);} /// dittopragma(inline, true)double tan(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) tan(cast(real) x) : tanImpl(x); } /// dittopragma(inline, true)float tan(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) tan(cast(real) x) : tanImpl(x); } ///@safe unittest{    import std.math.operations : isClose;    import std.math.traits : isIdentical;    import std.math.constants : PI;    import std.math.algebraic : sqrt;     assert(isIdentical(tan(0.0), 0.0));    assert(tan(PI).isClose(0, 0.0, 1e-10));    assert(tan(PI / 3).isClose(sqrt(3.0)));} version (InlineAsm_X87)private real tanAsm(real x) @trusted pure nothrow @nogc{    // Separating `return real.nan` from the asm block on LDC produces unintended    // behaviour as additional instructions are generated, invalidating the asm    // logic inside the previous block. To circumvent this, we can push rnan    // manually by creating an immutable variable in the stack.    immutable rnan = real.nan;     version (X86)    {    asm pure nothrow @nogc    {        fld     x[EBP]                  ; // load theta        fxam                            ; // test for oddball values        fstsw   AX                      ;        sahf                            ;        jc      trigerr                 ; // x is NAN, infinity, or empty                                          // 387's can handle subnormalsSC18:   fptan                           ;        fstsw   AX                      ;        sahf                            ;        jnp     Clear1                  ; // C2 = 1 (x is out of range)         // Do argument reduction to bring x into range        fldpi                           ;        fxch                            ;SC17:   fprem1                          ;        fstsw   AX                      ;        sahf                            ;        jp      SC17                    ;        fstp    ST(1)                   ; // remove pi from stack        jmp     SC18                    ; trigerr:        jnp     Lret                    ; // if theta is NAN, return theta        fstp    ST(0)                   ; // dump theta        fld     rnan                    ; // return rnan        jmp     Lret                    ;Clear1:        fstp    ST(0)                   ; // dump X, which is always 1Lret:        ;    }    }    else version (X86_64)    {        version (Win64)        {            asm pure nothrow @nogc            {                fld     real ptr [RCX]  ; // load theta            }        }        else        {            asm pure nothrow @nogc            {                fld     x[RBP]          ; // load theta            }        }    asm pure nothrow @nogc    {        fxam                            ; // test for oddball values        fstsw   AX                      ;        test    AH,1                    ;        jnz     trigerr                 ; // x is NAN, infinity, or empty                                          // 387's can handle subnormalsSC18:   fptan                           ;        fstsw   AX                      ;        test    AH,4                    ;        jz      Clear1                  ; // C2 = 1 (x is out of range)         // Do argument reduction to bring x into range        fldpi                           ;        fxch                            ;SC17:   fprem1                          ;        fstsw   AX                      ;        test    AH,4                    ;        jnz     SC17                    ;        fstp    ST(1)                   ; // remove pi from stack        jmp     SC18                    ; trigerr:        test    AH,4                    ;        jz      Lret                    ; // if theta is NAN, return theta        fstp    ST(0)                   ; // dump theta        fld     rnan                    ; // return rnan        jmp     Lret                    ;Clear1:        fstp    ST(0)                   ; // dump X, which is always 1Lret:        ;    }    }    else        static assert(0);} private T tanImpl(T)(T x) @safe pure nothrow @nogc{    import std.math.traits : floatTraits, RealFormat;    import std.math.constants : PI, PI_4;    import std.math.rounding : floor;    import std.math.algebraic : poly;    import std.math.traits : isInfinity, isNaN, signbit;     // Coefficients for tan(x) and PI/4 split into three parts.    enum realFormat = floatTraits!T.realFormat;    static if (realFormat == RealFormat.ieeeQuadruple)    {        static immutable T[6] P = [            2.883414728874239697964612246732416606301E10L,            -2.307030822693734879744223131873392503321E9L,            5.160188250214037865511600561074819366815E7L,            -4.249691853501233575668486667664718192660E5L,            1.272297782199996882828849455156962260810E3L,            -9.889929415807650724957118893791829849557E-1L        ];        static immutable T[7] Q = [            8.650244186622719093893836740197250197602E10L,            -4.152206921457208101480801635640958361612E10L,            2.758476078803232151774723646710890525496E9L,            -5.733709132766856723608447733926138506824E7L,            4.529422062441341616231663543669583527923E5L,            -1.317243702830553658702531997959756728291E3L,            1.0        ];         enum T P1 =            7.853981633974483067550664827649598009884357452392578125E-1L;        enum T P2 =            2.8605943630549158983813312792950660807511260829685741796657E-18L;        enum T P3 =            2.1679525325309452561992610065108379921905808E-35L;    }    else static if (realFormat == RealFormat.ieeeExtended ||                    realFormat == RealFormat.ieeeDouble)    {        static immutable T[3] P = [           -1.7956525197648487798769E7L,            1.1535166483858741613983E6L,           -1.3093693918138377764608E4L,        ];        static immutable T[5] Q = [           -5.3869575592945462988123E7L,            2.5008380182335791583922E7L,           -1.3208923444021096744731E6L,            1.3681296347069295467845E4L,            1.0000000000000000000000E0L,        ];         enum T P1 = 7.853981554508209228515625E-1L;        enum T P2 = 7.946627356147928367136046290398E-9L;        enum T P3 = 3.061616997868382943065164830688E-17L;    }    else static if (realFormat == RealFormat.ieeeSingle)    {        static immutable T[6] P = [            3.33331568548E-1,            1.33387994085E-1,            5.34112807005E-2,            2.44301354525E-2,            3.11992232697E-3,            9.38540185543E-3,        ];         enum T P1 = 0.78515625;        enum T P2 = 2.4187564849853515625E-4;        enum T P3 = 3.77489497744594108E-8;    }    else        static assert(0, "no coefficients for tan()");     // Special cases.    if (x == cast(T) 0.0 || isNaN(x))        return x;    if (isInfinity(x))        return T.nan;     // Make argument positive but save the sign.    bool sign = false;    if (signbit(x))    {        sign = true;        x = -x;    }     // Compute x mod PI/4.    static if (realFormat == RealFormat.ieeeSingle)    {        enum T FOPI = 4 / PI;        int j = cast(int) (FOPI * x);        T y = j;        T z;    }    else    {        T y = floor(x / cast(T) PI_4);        // Strip high bits of integer part.        enum T highBitsFactor = (realFormat == RealFormat.ieeeDouble ? 0x1p3 : 0x1p4);        enum T highBitsInv = 1.0 / highBitsFactor;        T z = y * highBitsInv;        // Compute y - 2^numHighBits * (y / 2^numHighBits).        z = y - highBitsFactor * floor(z);         // Integer and fraction part modulo one octant.        int j = cast(int)(z);    }     // Map zeros and singularities to origin.    if (j & 1)    {        j += 1;        y += cast(T) 1.0;    }     z = ((x - y * P1) - y * P2) - y * P3;    const T zz = z * z;     enum T zzThreshold = (realFormat == RealFormat.ieeeSingle ? 1.0e-4L :                          realFormat == RealFormat.ieeeDouble ? 1.0e-14L : 1.0e-20L);    if (zz > zzThreshold)    {        static if (realFormat == RealFormat.ieeeSingle)            y = z + z * (zz * poly(zz, P));        else            y = z + z * (zz * poly(zz, P) / poly(zz, Q));    }    else        y = z;     if (j & 2)        y = (cast(T) -1.0) / y;     return (sign) ? -y : y;} @safe @nogc nothrow unittest{    static void testTan(T)()    {        import std.math.operations : CommonDefaultFor, isClose, NaN;        import std.math.traits : isIdentical, isNaN;        import std.math.constants : PI, PI_4;         // ±0        const T zero = 0.0;        assert(isIdentical(tan(zero), zero));        assert(isIdentical(tan(-zero), -zero));        // ±∞        const T inf = T.infinity;        assert(isNaN(tan(inf)));        assert(isNaN(tan(-inf)));        // NaN        const T specialNaN = NaN(0x0123L);        assert(isIdentical(tan(specialNaN), specialNaN));         static immutable T[2][] vals =        [            // angle, tan            [   .5,  .546302489843790513255L],            [   1,   1.55740772465490223050L],            [   1.5, 14.1014199471717193876L],            [   2,  -2.18503986326151899164L],            [   2.5,-.747022297238660279355L],            [   3,  -.142546543074277805295L],            [   3.5, .374585640158594666330L],            [   4,   1.15782128234957758313L],            [   4.5, 4.63733205455118446831L],            [   5,  -3.38051500624658563698L],            [   5.5,-.995584052213885017701L],            [   6,  -.291006191384749157053L],            [   6.5, .220277200345896811825L],            [   10,  .648360827459086671259L],             // special angles            [   PI_4,   1],            //[   PI_2,   T.infinity], // PI_2 is not _exactly_ pi/2.            [   3*PI_4, -1],            [   PI,     0],            [   5*PI_4, 1],            //[   3*PI_2, -T.infinity],            [   7*PI_4, -1],            [   2*PI,   0],         ];         foreach (ref val; vals)        {            T x = val[0];            T r = val[1];            T t = tan(x);

Lines 1–500 of 1424.

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