math.d - This D code implements mathematical functions and …

/std/math.d

http://github.com/jcd/phobos · D · 5827 lines · 3914 code · 563 blank · 1350 comment · 779 complexity · 9fbdeb980627c624e3f89cb331d11a13 MD5 · raw file
Large files are truncated click here to view the full file

// Written in the D programming language.

/**
 * Elementary mathematical functions
 *
 * Contains the elementary mathematical functions (powers, roots,
 * and trignometric functions), and low-level floating-point operations.
 * Mathematical special functions are available in std.mathspecial.
 *
 * The functionality closely follows the IEEE754-2008 standard for
 * floating-point arithmetic, including the use of camelCase names rather
 * than C99-style lower case names. All of these functions behave correctly
 * when presented with an infinity or NaN.
 *
 * Unlike C, there is no global 'errno' variable. Consequently, almost all of
 * these functions are pure nothrow.
 *
 * Status:
 * The semantics and names of feqrel and approxEqual will be revised.
 *
 * Macros:
 *      WIKI = Phobos/StdMath
 *
 *      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)
 *      SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
 *      GAMMA = &#915;
 *      THETA = &theta;
 *      INTEGRAL = &#8747;
 *      INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
 *      POWER = $1<sup>$2</sup>
 *      SUB = $1<sub>$2</sub>
 *      BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
 *      CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
 *      PLUSMN = &plusmn;
 *      INFIN = &infin;
 *      PLUSMNINF = &plusmn;&infin;
 *      PI = &pi;
 *      LT = &lt;
 *      GT = &gt;
 *      SQRT = &radic;
 *      HALF = &frac12;
 *
 * Copyright: Copyright Digital Mars 2000 - 2011.
 *            D implementations of tan, atan, atan2, exp, expm1, exp2, log, log10, log1p,
 *            log2, floor, ceil and lrint functions are based on the CEPHES math library,
 *            which is Copyright (C) 2001 Stephen L. Moshier <steve@moshier.net>
 *            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:   <a href="http://www.boost.org/LICENSE_1_0.txt">Boost License 1.0</a>.
 * Authors:   $(WEB digitalmars.com, Walter Bright),
 *                        Don Clugston, Conversion of CEPHES math library to D by Iain Buclaw
 * Source: $(PHOBOSSRC std/_math.d)
 */
module std.math;

import core.stdc.math;
import std.traits;

version(unittest)
{
    import std.typetuple;
}

version(LDC)
{
    import ldc.intrinsics;
}

version(DigitalMars)
{
    version = INLINE_YL2X;        // x87 has opcodes for these
}

version (X86)
{
    version = X86_Any;
}

version (X86_64)
{
    version = X86_Any;
}

version(D_InlineAsm_X86)
{
    version = InlineAsm_X86_Any;
}
else version(D_InlineAsm_X86_64)
{
    version = InlineAsm_X86_Any;
}


version(unittest)
{
    import core.stdc.stdio;

    static if(real.sizeof > double.sizeof)
        enum uint useDigits = 16;
    else
        enum uint useDigits = 15;

    /******************************************
     * Compare floating point numbers to n decimal digits of precision.
     * Returns:
     *  1       match
     *  0       nomatch
     */

    private bool equalsDigit(real x, real y, uint ndigits)
    {
        if (signbit(x) != signbit(y))
            return 0;

        if (isinf(x) && isinf(y))
            return 1;
        if (isinf(x) || isinf(y))
            return 0;

        if (isnan(x) && isnan(y))
            return 1;
        if (isnan(x) || isnan(y))
            return 0;

        char[30] bufx;
        char[30] bufy;
        assert(ndigits < bufx.length);

        int ix;
        int iy;
        ix = sprintf(bufx.ptr, "%.*Lg", ndigits, x);
        assert(ix < bufx.length && ix > 0);
        iy = sprintf(bufy.ptr, "%.*Lg", ndigits, y);
        assert(ix < bufy.length && ix > 0);

        return bufx[0 .. ix] == bufy[0 .. iy];
    }
}



private:
/*
 * The following IEEE 'real' formats are currently supported:
 * 64 bit Big-endian  'double' (eg PowerPC)
 * 128 bit Big-endian 'quadruple' (eg SPARC)
 * 64 bit Little-endian 'double' (eg x86-SSE2)
 * 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium).
 * 128 bit Little-endian 'quadruple' (not implemented on any known processor!)
 *
 * Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support
 */
version(LittleEndian)
{
    static assert(real.mant_dig == 53 || real.mant_dig==64
               || real.mant_dig == 113,
      "Only 64-bit, 80-bit, and 128-bit reals"
      " are supported for LittleEndian CPUs");
}
else
{
    static assert(real.mant_dig == 53 || real.mant_dig==106
               || real.mant_dig == 113,
    "Only 64-bit and 128-bit reals are supported for BigEndian CPUs."
    " double-double reals have partial support");
}

// Constants used for extracting the components of the representation.
// They supplement the built-in floating point properties.
template floatTraits(T)
{
    // EXPMASK is a ushort mask to select the exponent portion (without sign)
    // 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 T RECIP_EPSILON = (1/T.epsilon);
    static if (T.mant_dig == 24)
    { // float
        enum ushort EXPMASK = 0x7F80;
        enum ushort EXPBIAS = 0x3F00;
        enum uint EXPMASK_INT = 0x7F80_0000;
        enum uint MANTISSAMASK_INT = 0x007F_FFFF;
        version(LittleEndian)
        {
            enum EXPPOS_SHORT = 1;
        }
        else
        {
            enum EXPPOS_SHORT = 0;
        }
    }
    else static if (T.mant_dig == 53) // double, or real==double
    {
        enum ushort EXPMASK = 0x7FF0;
        enum ushort EXPBIAS = 0x3FE0;
        enum uint EXPMASK_INT = 0x7FF0_0000;
        enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only
        version(LittleEndian)
        {
            enum EXPPOS_SHORT = 3;
            enum SIGNPOS_BYTE = 7;
        }
        else
        {
            enum EXPPOS_SHORT = 0;
            enum SIGNPOS_BYTE = 0;
        }
    }
    else static if (T.mant_dig == 64) // real80
    {
        enum ushort EXPMASK = 0x7FFF;
        enum ushort EXPBIAS = 0x3FFE;
        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
    {
        enum ushort EXPMASK = 0x7FFF;
        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) // doubledouble
    {
        enum ushort EXPMASK = 0x7FF0;
        // the exponent byte is not unique
        version(LittleEndian)
        {
            enum EXPPOS_SHORT = 7; // [3] is also an exp short
            enum SIGNPOS_BYTE = 15;
        }
        else
        {
            enum EXPPOS_SHORT = 0; // [4] is also an exp short
            enum SIGNPOS_BYTE = 0;
        }
    }
}

// These apply to all floating-point types
version(LittleEndian)
{
    enum MANTISSA_LSB = 0;
    enum MANTISSA_MSB = 1;
}
else
{
    enum MANTISSA_LSB = 1;
    enum MANTISSA_MSB = 0;
}

public:

// Values obtained from Wolfram Alpha. 116 bits ought to be enough for anybody.
// Wolfram Alpha LLC. 2011. Wolfram|Alpha. http://www.wolframalpha.com/input/?i=e+in+base+16 (access July 6, 2011).
enum real E =          0x1.5bf0a8b1457695355fb8ac404e7a8p+1L; /** e = 2.718281... */
enum real LOG2T =      0x1.a934f0979a3715fc9257edfe9b5fbp+1L; /** $(SUB log, 2)10 = 3.321928... */
enum real LOG2E =      0x1.71547652b82fe1777d0ffda0d23a8p+0L; /** $(SUB log, 2)e = 1.442695... */
enum real LOG2 =       0x1.34413509f79fef311f12b35816f92p-2L; /** $(SUB log, 10)2 = 0.301029... */
enum real LOG10E =     0x1.bcb7b1526e50e32a6ab7555f5a67cp-2L; /** $(SUB log, 10)e = 0.434294... */
enum real LN2 =        0x1.62e42fefa39ef35793c7673007e5fp-1L; /** ln 2  = 0.693147... */
enum real LN10 =       0x1.26bb1bbb5551582dd4adac5705a61p+1L; /** ln 10 = 2.302585... */
enum real PI =         0x1.921fb54442d18469898cc51701b84p+1L; /** $(_PI) = 3.141592... */
enum real PI_2 =       PI/2;                                  /** $(PI) / 2 = 1.570796... */
enum real PI_4 =       PI/4;                                  /** $(PI) / 4 = 0.785398... */
enum real M_1_PI =     0x1.45f306dc9c882a53f84eafa3ea69cp-2L; /** 1 / $(PI) = 0.318309... */
enum real M_2_PI =     2*M_1_PI;                              /** 2 / $(PI) = 0.636619... */
enum real M_2_SQRTPI = 0x1.20dd750429b6d11ae3a914fed7fd8p+0L; /** 2 / $(SQRT)$(PI) = 1.128379... */
enum real SQRT2 =      0x1.6a09e667f3bcc908b2fb1366ea958p+0L; /** $(SQRT)2 = 1.414213... */
enum real SQRT1_2 =    SQRT2/2;                               /** $(SQRT)$(HALF) = 0.707106... */
// Note: Make sure the magic numbers in compiler backend for x87 match these.


/***********************************
 * Calculates the absolute value
 *
 * For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) )
 * = hypot(z.re, z.im).
 */
Num abs(Num)(Num x) @safe pure nothrow
    if (is(typeof(Num.init >= 0)) && is(typeof(-Num.init)) &&
            !(is(Num* : const(ifloat*)) || is(Num* : const(idouble*))
                    || is(Num* : const(ireal*))))
{
    static if (isFloatingPoint!(Num))
        return fabs(x);
    else
        return x>=0 ? x : -x;
}

auto abs(Num)(Num z) @safe pure nothrow
    if (is(Num* : const(cfloat*)) || is(Num* : const(cdouble*))
            || is(Num* : const(creal*)))
{
    return hypot(z.re, z.im);
}

/** ditto */
real abs(Num)(Num y) @safe pure nothrow
    if (is(Num* : const(ifloat*)) || is(Num* : const(idouble*))
            || is(Num* : const(ireal*)))
{
    return fabs(y.im);
}


unittest
{
    assert(isIdentical(abs(-0.0L), 0.0L));
    assert(isNaN(abs(real.nan)));
    assert(abs(-real.infinity) == real.infinity);
    assert(abs(-3.2Li) == 3.2L);
    assert(abs(71.6Li) == 71.6L);
    assert(abs(-56) == 56);
    assert(abs(2321312L)  == 2321312L);
    assert(abs(-1+1i) == sqrt(2.0L));
}

/***********************************
 * Complex conjugate
 *
 *  conj(x + iy) = x - iy
 *
 * Note that z * conj(z) = $(POWER z.re, 2) - $(POWER z.im, 2)
 * is always a real number
 */
creal conj(creal z) @safe pure nothrow
{
    return z.re - z.im*1i;
}

/** ditto */
ireal conj(ireal y) @safe pure nothrow
{
    return -y;
}


unittest
{
    assert(conj(7 + 3i) == 7-3i);
    ireal z = -3.2Li;
    assert(conj(z) == -z);
}

/***********************************
 * 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).
 */

real cos(real x) @safe pure nothrow;       /* intrinsic */

/***********************************
 * Returns sine of x. x is in radians.
 *
 *      $(TABLE_SV
 *      $(TR $(TH x)               $(TH sin(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))
 *      )
 * Bugs:
 *      Results are undefined if |x| >= $(POWER 2,64).
 */

real sin(real x) @safe pure nothrow;       /* intrinsic */


/***********************************
 *  sine, complex and imaginary
 *
 *  sin(z) = sin(z.re)*cosh(z.im) + cos(z.re)*sinh(z.im)i
 *
 * If both sin($(THETA)) and cos($(THETA)) are required,
 * it is most efficient to use expi($(THETA)).
 */
creal sin(creal z) @safe pure nothrow
{
    creal cs = expi(z.re);
    creal csh = coshisinh(z.im);
    return cs.im * csh.re + cs.re * csh.im * 1i;
}

/** ditto */
ireal sin(ireal y) @safe pure nothrow
{
    return cosh(y.im)*1i;
}

unittest
{
  assert(sin(0.0+0.0i) == 0.0);
  assert(sin(2.0+0.0i) == sin(2.0L) );
}

/***********************************
 *  cosine, complex and imaginary
 *
 *  cos(z) = cos(z.re)*cosh(z.im) - sin(z.re)*sinh(z.im)i
 */
creal cos(creal z) @safe pure nothrow
{
    creal cs = expi(z.re);
    creal csh = coshisinh(z.im);
    return cs.re * csh.re - cs.im * csh.im * 1i;
}

/** ditto */
real cos(ireal y) @safe pure nothrow
{
    return cosh(y.im);
}

unittest
{
    assert(cos(0.0+0.0i)==1.0);
    assert(cos(1.3L+0.0i)==cos(1.3L));
    assert(cos(5.2Li)== cosh(5.2L));
}

/****************************************************************************
 * 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))
 *      )
 */

real tan(real x) @trusted pure nothrow
{
    version(D_InlineAsm_X86)
    {
    asm
    {
        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 subnormals
SC18:   fptan                           ;
        fstp    ST(0)                   ; // dump X, which is always 1
        fstsw   AX                      ;
        sahf                            ;
        jnp     Lret                    ; // 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
    }
    return real.nan;

Lret: {}
    }
    else version(D_InlineAsm_X86_64)
    {
        version (Win64)
        {
            asm
            {
                fld     real ptr [RCX]  ; // load theta
            }
        }
        else
        {
            asm
            {
                fld     x[RBP]          ; // load theta
            }
        }
    asm
    {
        fxam                            ; // test for oddball values
        fstsw   AX                      ;
        test    AH,1                    ;
        jnz     trigerr                 ; // x is NAN, infinity, or empty
                                          // 387's can handle subnormals
SC18:   fptan                           ;
        fstp    ST(0)                   ; // dump X, which is always 1
        fstsw   AX                      ;
        test    AH,4                    ;
        jz      Lret                    ; // 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
    }
    return real.nan;

Lret: {}
    }
    else
    {
        // Coefficients for tan(x)
        static immutable real[3] P = [
           -1.7956525197648487798769E7L,
            1.1535166483858741613983E6L,
           -1.3093693918138377764608E4L,
        ];
        static immutable real[5] Q = [
           -5.3869575592945462988123E7L,
            2.5008380182335791583922E7L,
           -1.3208923444021096744731E6L,
            1.3681296347069295467845E4L,
            1.0000000000000000000000E0L,
        ];

        // PI/4 split into three parts.
        enum real P1 = 7.853981554508209228515625E-1L;
        enum real P2 = 7.946627356147928367136046290398E-9L;
        enum real P3 = 3.061616997868382943065164830688E-17L;

        // Special cases.
        if (x == 0.0 || isNaN(x))
            return x;
        if (isInfinity(x))
            return real.nan;

        // Make argument positive but save the sign.
        bool sign = false;
        if (signbit(x))
        {
            sign = true;
            x = -x;
        }

        // Compute x mod PI/4.
        real y = floor(x / PI_4);
        // Strip high bits of integer part.
        real z = ldexp(y, -4);
        // Compute y - 16 * (y / 16).
        z = y - ldexp(floor(z), 4);

        // Integer and fraction part modulo one octant.
        int j = cast(int)(z);

        // Map zeros and singularities to origin.
        if (j & 1)
        {
            j += 1;
            y += 1.0;
        }

        z = ((x - y * P1) - y * P2) - y * P3;
        real zz = z * z;

        if (zz > 1.0e-20L)
            y = z + z * (zz * poly(zz, P) / poly(zz, Q));
        else
            y = z;

        if (j & 2)
            y = -1.0 / y;

        return (sign) ? -y : y;
    }
}

unittest
{
    static real[2][] vals =     // angle,tan
        [
         [   0,   0],
         [   .5,  .5463024898],
         [   1,   1.557407725],
         [   1.5, 14.10141995],
         [   2,  -2.185039863],
         [   2.5,-.7470222972],
         [   3,  -.1425465431],
         [   3.5, .3745856402],
         [   4,   1.157821282],
         [   4.5, 4.637332055],
         [   5,  -3.380515006],
         [   5.5,-.9955840522],
         [   6,  -.2910061914],
         [   6.5, .2202772003],
         [   10,  .6483608275],

         // special angles
         [   PI_4,   1],
         //[   PI_2,   real.infinity], // PI_2 is not _exactly_ pi/2.
         [   3*PI_4, -1],
         [   PI,     0],
         [   5*PI_4, 1],
         //[   3*PI_2, -real.infinity],
         [   7*PI_4, -1],
         [   2*PI,   0],
         ];
    int i;

    for (i = 0; i < vals.length; i++)
    {
        real x = vals[i][0];
        real r = vals[i][1];
        real t = tan(x);

        //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
        if (!isIdentical(r, t)) assert(fabs(r-t) <= .0000001);

        x = -x;
        r = -r;
        t = tan(x);
        //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
        if (!isIdentical(r, t) && !(r!=r && t!=t)) assert(fabs(r-t) <= .0000001);
    }
    // overflow
    assert(isNaN(tan(real.infinity)));
    assert(isNaN(tan(-real.infinity)));
    // NaN propagation
    assert(isIdentical( tan(NaN(0x0123L)), NaN(0x0123L) ));
}

unittest
{
    assert(equalsDigit(tan(PI / 3), std.math.sqrt(3.0), useDigits));
}

/***************
 * Calculates the arc cosine of x,
 * returning a value ranging from 0 to $(PI).
 *
 *      $(TABLE_SV
 *      $(TR $(TH x)         $(TH acos(x)) $(TH invalid?))
 *      $(TR $(TD $(GT)1.0)  $(TD $(NAN))  $(TD yes))
 *      $(TR $(TD $(LT)-1.0) $(TD $(NAN))  $(TD yes))
 *      $(TR $(TD $(NAN))    $(TD $(NAN))  $(TD yes))
 *  )
 */
real acos(real x) @safe pure nothrow
{
    return atan2(sqrt(1-x*x), x);
}

/// ditto
double acos(double x) @safe pure nothrow { return acos(cast(real)x); }

/// ditto
float acos(float x) @safe pure nothrow  { return acos(cast(real)x); }

unittest
{
    assert(equalsDigit(acos(0.5), std.math.PI / 3, useDigits));
}

/***************
 * Calculates the arc sine of x,
 * returning a value ranging from -$(PI)/2 to $(PI)/2.
 *
 *      $(TABLE_SV
 *      $(TR $(TH x)            $(TH asin(x))      $(TH invalid?))
 *      $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
 *      $(TR $(TD $(GT)1.0)     $(TD $(NAN))       $(TD yes))
 *      $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD yes))
 *  )
 */
real asin(real x) @safe pure nothrow
{
    return atan2(x, sqrt(1-x*x));
}

/// ditto
double asin(double x) @safe pure nothrow { return asin(cast(real)x); }

/// ditto
float asin(float x) @safe pure nothrow  { return asin(cast(real)x); }

unittest
{
    assert(equalsDigit(asin(0.5), PI / 6, useDigits));
}

/***************
 * Calculates the arc tangent of x,
 * returning a value ranging from -$(PI)/2 to $(PI)/2.
 *
 *      $(TABLE_SV
 *      $(TR $(TH x)                 $(TH atan(x))      $(TH invalid?))
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no))
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))       $(TD yes))
 *  )
 */
real atan(real x) @safe pure nothrow
{
    version(InlineAsm_X86_Any)
    {
        return atan2(x, 1.0L);
    }
    else
    {
        // Coefficients for atan(x)
        static immutable real[5] P = [
           -5.0894116899623603312185E1L,
           -9.9988763777265819915721E1L,
           -6.3976888655834347413154E1L,
           -1.4683508633175792446076E1L,
           -8.6863818178092187535440E-1L,
        ];
        static immutable real[6] Q = [
            1.5268235069887081006606E2L,
            3.9157570175111990631099E2L,
            3.6144079386152023162701E2L,
            1.4399096122250781605352E2L,
            2.2981886733594175366172E1L,
            1.0000000000000000000000E0L,
        ];

        // tan(PI/8)
        enum real TAN_PI_8 = 4.1421356237309504880169e-1L;
        // tan(3 * PI/8)
        enum real TAN3_PI_8 = 2.41421356237309504880169L;

        // Special cases.
        if (x == 0.0)
            return x;
        if (isInfinity(x))
            return copysign(PI_2, x);

        // Make argument positive but save the sign.
        bool sign = false;
        if (signbit(x))
        {
            sign = true;
            x = -x;
        }

        // Range reduction.
        real y;
        if (x > TAN3_PI_8)
        {
            y = PI_2;
            x = -(1.0 / x);
        }
        else if (x > TAN_PI_8)
        {
            y = PI_4;
            x = (x - 1.0)/(x + 1.0);
        }
        else
            y = 0.0;

        // Rational form in x^^2.
        real z = x * x;
        y = y + (poly(z, P) / poly(z, Q)) * z * x + x;

        return (sign) ? -y : y;
    }
}

/// ditto
double atan(double x) @safe pure nothrow { return atan(cast(real)x); }

/// ditto
float atan(float x)  @safe pure nothrow { return atan(cast(real)x); }

unittest
{
    assert(equalsDigit(atan(std.math.sqrt(3.0)), PI / 3, useDigits));
}

/***************
 * Calculates the arc tangent of y / x,
 * returning a value ranging from -$(PI) to $(PI).
 *
 *      $(TABLE_SV
 *      $(TR $(TH y)                 $(TH x)            $(TH atan(y, x)))
 *      $(TR $(TD $(NAN))            $(TD anything)     $(TD $(NAN)) )
 *      $(TR $(TD anything)          $(TD $(NAN))       $(TD $(NAN)) )
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(GT)0.0)     $(TD $(PLUSMN)0.0) )
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0)         $(TD $(PLUSMN)0.0) )
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(LT)0.0)     $(TD $(PLUSMN)$(PI)))
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD -0.0)         $(TD $(PLUSMN)$(PI)))
 *      $(TR $(TD $(GT)0.0)          $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
 *      $(TR $(TD $(LT)0.0)          $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) )
 *      $(TR $(TD $(GT)0.0)          $(TD $(INFIN))     $(TD $(PLUSMN)0.0) )
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything)     $(TD $(PLUSMN)$(PI)/2))
 *      $(TR $(TD $(GT)0.0)          $(TD -$(INFIN))    $(TD $(PLUSMN)$(PI)) )
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN))     $(TD $(PLUSMN)$(PI)/4))
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN))    $(TD $(PLUSMN)3$(PI)/4))
 *      )
 */
real atan2(real y, real x) @trusted pure nothrow
{
    version(InlineAsm_X86_Any)
    {
        version (Win64)
        {
            asm {
                naked;
                fld real ptr [RDX]; // y
                fld real ptr [RCX]; // x
                fpatan;
                ret;
            }
        }
        else
        {
            asm {
                fld y;
                fld x;
                fpatan;
            }
        }
    }
    else
    {
        // Special cases.
        if (isNaN(x) || isNaN(y))
            return real.nan;
        if (y == 0.0)
        {
            if (x >= 0 && !signbit(x))
                return copysign(0, y);
            else
                return copysign(PI, y);
        }
        if (x == 0.0)
            return copysign(PI_2, y);
        if (isInfinity(x))
        {
            if (signbit(x))
            {
                if (isInfinity(y))
                    return copysign(3*PI_4, y);
                else
                    return copysign(PI, y);
            }
            else
            {
                if (isInfinity(y))
                    return copysign(PI_4, y);
                else
                    return copysign(0.0, y);
            }
        }
        if (isInfinity(y))
            return copysign(PI_2, y);

        // Call atan and determine the quadrant.
        real z = atan(y / x);

        if (signbit(x))
        {
            if (signbit(y))
                z = z - PI;
            else
                z = z + PI;
        }

        if (z == 0.0)
            return copysign(z, y);

        return z;
    }
}

/// ditto
double atan2(double y, double x) @safe pure nothrow
{
    return atan2(cast(real)y, cast(real)x);
}

/// ditto
float atan2(float y, float x) @safe pure nothrow
{
    return atan2(cast(real)y, cast(real)x);
}

unittest
{
    assert(equalsDigit(atan2(1.0L, std.math.sqrt(3.0L)), PI / 6, useDigits));
}

/***********************************
 * Calculates the hyperbolic cosine of x.
 *
 *      $(TABLE_SV
 *      $(TR $(TH x)                 $(TH cosh(x))      $(TH invalid?))
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
 *      )
 */
real cosh(real x) @safe pure nothrow
{
    //  cosh = (exp(x)+exp(-x))/2.
    // The naive implementation works correctly.
    real y = exp(x);
    return (y + 1.0/y) * 0.5;
}

/// ditto
double cosh(double x) @safe pure nothrow { return cosh(cast(real)x); }

/// ditto
float cosh(float x) @safe pure nothrow  { return cosh(cast(real)x); }

unittest
{
    assert(equalsDigit(cosh(1.0), (E + 1.0 / E) / 2, useDigits));
}

/***********************************
 * Calculates the hyperbolic sine of x.
 *
 *      $(TABLE_SV
 *      $(TR $(TH x)                 $(TH sinh(x))           $(TH invalid?))
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no))
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
 *      )
 */
real sinh(real x) @safe pure nothrow
{
    //  sinh(x) =  (exp(x)-exp(-x))/2;
    // Very large arguments could cause an overflow, but
    // the maximum value of x for which exp(x) + exp(-x)) != exp(x)
    // is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80.
    if (fabs(x) > real.mant_dig * LN2)
    {
        return copysign(0.5 * exp(fabs(x)), x);
    }

    real y = expm1(x);
    return 0.5 * y / (y+1) * (y+2);
}

/// ditto
double sinh(double x) @safe pure nothrow { return sinh(cast(real)x); }

/// ditto
float sinh(float x) @safe pure nothrow  { return sinh(cast(real)x); }

unittest
{
    assert(equalsDigit(sinh(1.0), (E - 1.0 / E) / 2, useDigits));
}

/***********************************
 * Calculates the hyperbolic tangent of x.
 *
 *      $(TABLE_SV
 *      $(TR $(TH x)                 $(TH tanh(x))      $(TH invalid?))
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no) )
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
 *      )
 */
real tanh(real x) @safe pure nothrow
{
    //  tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
    if (fabs(x) > real.mant_dig * LN2)
    {
        return copysign(1, x);
    }

    real y = expm1(2*x);
    return y / (y + 2);
}

/// ditto
double tanh(double x) @safe pure nothrow { return tanh(cast(real)x); }

/// ditto
float tanh(float x) @safe pure nothrow { return tanh(cast(real)x); }

unittest
{
    assert(equalsDigit(tanh(1.0), sinh(1.0) / cosh(1.0), 15));
}

package:

/* Returns cosh(x) + I * sinh(x)
 * Only one call to exp() is performed.
 */
creal coshisinh(real x) @safe pure nothrow
{
    // See comments for cosh, sinh.
    if (fabs(x) > real.mant_dig * LN2)
    {
        real y = exp(fabs(x));
        return y * 0.5 + 0.5i * copysign(y, x);
    }
    else
    {
        real y = expm1(x);
        return (y + 1.0 + 1.0/(y + 1.0)) * 0.5 + 0.5i * y / (y+1) * (y+2);
    }
}

unittest
{
    creal c = coshisinh(3.0L);
    assert(c.re == cosh(3.0L));
    assert(c.im == sinh(3.0L));
}

public:

/***********************************
 * Calculates the inverse hyperbolic cosine of x.
 *
 *  Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
 *
 * $(TABLE_DOMRG
 *  $(DOMAIN 1..$(INFIN))
 *  $(RANGE  1..log(real.max), $(INFIN)) )
 *      $(TABLE_SV
 *    $(SVH  x,     acosh(x) )
 *    $(SV  $(NAN), $(NAN) )
 *    $(SV  $(LT)1,     $(NAN) )
 *    $(SV  1,      0       )
 *    $(SV  +$(INFIN),+$(INFIN))
 *  )
 */
real acosh(real x) @safe pure nothrow
{
    if (x > 1/real.epsilon)
        return LN2 + log(x);
    else
        return log(x + sqrt(x*x - 1));
}

/// ditto
double acosh(double x) @safe pure nothrow { return acosh(cast(real)x); }

/// ditto
float acosh(float x) @safe pure nothrow  { return acosh(cast(real)x); }


unittest
{
    assert(isNaN(acosh(0.9)));
    assert(isNaN(acosh(real.nan)));
    assert(acosh(1.0)==0.0);
    assert(acosh(real.infinity) == real.infinity);
    assert(isNaN(acosh(0.5)));
    assert(equalsDigit(acosh(cosh(3.0)), 3, useDigits));
}

/***********************************
 * Calculates the inverse hyperbolic sine of x.
 *
 *  Mathematically,
 *  ---------------
 *  asinh(x) =  log( x + sqrt( x*x + 1 )) // if x >= +0
 *  asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
 *  -------------
 *
 *    $(TABLE_SV
 *    $(SVH x,                asinh(x)       )
 *    $(SV  $(NAN),           $(NAN)         )
 *    $(SV  $(PLUSMN)0,       $(PLUSMN)0      )
 *    $(SV  $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
 *    )
 */
real asinh(real x) @safe pure nothrow
{
    return (fabs(x) > 1 / real.epsilon)
       // beyond this point, x*x + 1 == x*x
       ?  copysign(LN2 + log(fabs(x)), x)
       // sqrt(x*x + 1) ==  1 + x * x / ( 1 + sqrt(x*x + 1) )
       : copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
}

/// ditto
double asinh(double x) @safe pure nothrow { return asinh(cast(real)x); }

/// ditto
float asinh(float x) @safe pure nothrow { return asinh(cast(real)x); }

unittest
{
    assert(isIdentical(asinh(0.0), 0.0));
    assert(isIdentical(asinh(-0.0), -0.0));
    assert(asinh(real.infinity) == real.infinity);
    assert(asinh(-real.infinity) == -real.infinity);
    assert(isNaN(asinh(real.nan)));
    assert(equalsDigit(asinh(sinh(3.0)), 3, useDigits));
}

/***********************************
 * Calculates the inverse hyperbolic tangent of x,
 * returning a value from ranging from -1 to 1.
 *
 * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
 *
 *
 * $(TABLE_DOMRG
 *  $(DOMAIN -$(INFIN)..$(INFIN))
 *  $(RANGE  -1..1) )
 * $(TABLE_SV
 *    $(SVH  x,     acosh(x) )
 *    $(SV  $(NAN), $(NAN) )
 *    $(SV  $(PLUSMN)0, $(PLUSMN)0)
 *    $(SV  -$(INFIN), -0)
 * )
 */
real atanh(real x) @safe pure nothrow
{
    // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
    return  0.5 * log1p( 2 * x / (1 - x) );
}

/// ditto
double atanh(double x) @safe pure nothrow { return atanh(cast(real)x); }

/// ditto
float atanh(float x) @safe pure nothrow { return atanh(cast(real)x); }


unittest
{
    assert(isIdentical(atanh(0.0), 0.0));
    assert(isIdentical(atanh(-0.0),-0.0));
    assert(isNaN(atanh(real.nan)));
    assert(isNaN(atanh(-real.infinity)));
    assert(atanh(0.0) == 0);
    assert(equalsDigit(atanh(tanh(0.5L)), 0.5, useDigits));
}

/*****************************************
 * Returns x rounded to a long value using the current rounding mode.
 * If the integer value of x is
 * greater than long.max, the result is
 * indeterminate.
 */
long rndtol(real x) @safe pure nothrow;    /* intrinsic */


/*****************************************
 * Returns x rounded to a long value using the FE_TONEAREST rounding mode.
 * If the integer value of x is
 * greater than long.max, the result is
 * indeterminate.
 */
extern (C) real rndtonl(real x);

/***************************************
 * Compute square root of x.
 *
 *      $(TABLE_SV
 *      $(TR $(TH x)         $(TH sqrt(x))   $(TH invalid?))
 *      $(TR $(TD -0.0)      $(TD -0.0)      $(TD no))
 *      $(TR $(TD $(LT)0.0)  $(TD $(NAN))    $(TD yes))
 *      $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
 *      )
 */
float sqrt(float x) @safe pure nothrow;    /* intrinsic */

/// ditto
double sqrt(double x) @safe pure nothrow;  /* intrinsic */

/// ditto
real sqrt(real x) @safe pure nothrow;      /* intrinsic */

unittest
{
    //ctfe
    enum ZX80 = sqrt(7.0f);
    enum ZX81 = sqrt(7.0);
    enum ZX82 = sqrt(7.0L);
}

creal sqrt(creal z) @safe pure nothrow
{
    creal c;
    real x,y,w,r;

    if (z == 0)
    {
        c = 0 + 0i;
    }
    else
    {
        real z_re = z.re;
        real z_im = z.im;

        x = fabs(z_re);
        y = fabs(z_im);
        if (x >= y)
        {
            r = y / x;
            w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
        }
        else
        {
            r = x / y;
            w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
        }

        if (z_re >= 0)
        {
            c = w + (z_im / (w + w)) * 1.0i;
        }
        else
        {
            if (z_im < 0)
                w = -w;
            c = z_im / (w + w) + w * 1.0i;
        }
    }
    return c;
}

/**
 * Calculates e$(SUP x).
 *
 *  $(TABLE_SV
 *    $(TR $(TH x)             $(TH e$(SUP x)) )
 *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
 *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
 *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
 *  )
 */
real exp(real x) @trusted pure nothrow
{
    version(D_InlineAsm_X86)
    {
        //  e^^x = 2^^(LOG2E*x)
        // (This is valid because the overflow & underflow limits for exp
        // and exp2 are so similar).
        return exp2(LOG2E*x);
    }
    else version(D_InlineAsm_X86_64)
    {
        //  e^^x = 2^^(LOG2E*x)
        // (This is valid because the overflow & underflow limits for exp
        // and exp2 are so similar).
        return exp2(LOG2E*x);
    }
    else
    {
        // Coefficients for exp(x)
        static immutable real[3] P = [
            9.9999999999999999991025E-1L,
            3.0299440770744196129956E-2L,
            1.2617719307481059087798E-4L,
        ];
        static immutable real[4] Q = [
            2.0000000000000000000897E0L,
            2.2726554820815502876593E-1L,
            2.5244834034968410419224E-3L,
            3.0019850513866445504159E-6L,
        ];

        // C1 + C2 = LN2.
        enum real C1 = 6.9314575195312500000000E-1L;
        enum real C2 = 1.428606820309417232121458176568075500134E-6L;

        // Overflow and Underflow limits.
        enum real OF =  1.1356523406294143949492E4L;
        enum real UF = -1.1432769596155737933527E4L;

        // Special cases.
        if (isNaN(x))
            return x;
        if (x > OF)
            return real.infinity;
        if (x < UF)
            return 0.0;

        // Express: e^^x = e^^g * 2^^n
        //   = e^^g * e^^(n * LOG2E)
        //   = e^^(g + n * LOG2E)
        int n = cast(int)floor(LOG2E * x + 0.5);
        x -= n * C1;
        x -= n * C2;

        // Rational approximation for exponential of the fractional part:
        //  e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
        real xx = x * x;
        real px = x * poly(xx, P);
        x = px / (poly(xx, Q) - px);
        x = 1.0 + ldexp(x, 1);

        // Scale by power of 2.
        x = ldexp(x, n);

        return x;
    }
}

/// ditto
double exp(double x) @safe pure nothrow  { return exp(cast(real)x); }

/// ditto
float exp(float x)  @safe pure nothrow   { return exp(cast(real)x); }

unittest
{
    assert(equalsDigit(exp(3.0L), E * E * E, useDigits));
}

/**
 * Calculates the value of the natural logarithm base (e)
 * raised to the power of x, minus 1.
 *
 * For very small x, expm1(x) is more accurate
 * than exp(x)-1.
 *
 *  $(TABLE_SV
 *    $(TR $(TH x)             $(TH e$(SUP x)-1)  )
 *    $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) )
 *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN))    )
 *    $(TR $(TD -$(INFIN))     $(TD -1.0)         )
 *    $(TR $(TD $(NAN))        $(TD $(NAN))       )
 *  )
 */
real expm1(real x) @trusted pure nothrow
{
    version(D_InlineAsm_X86)
    {
        enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
        asm
        {
            /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
             * Author: Don Clugston.
             *
             *    expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
             *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
             *     and 2ym1 = (2^^(y-rndint(y))-1).
             *    If 2rndy  < 0.5*real.epsilon, result is -1.
             *    Implementation is otherwise the same as for exp2()
             */
            naked;
            fld real ptr [ESP+4] ; // x
            mov AX, [ESP+4+8]; // AX = exponent and sign
            sub ESP, 12+8; // Create scratch space on the stack
            // [ESP,ESP+2] = scratchint
            // [ESP+4..+6, +8..+10, +10] = scratchreal
            // set scratchreal mantissa = 1.0
            mov dword ptr [ESP+8], 0;
            mov dword ptr [ESP+8+4], 0x80000000;
            and AX, 0x7FFF; // drop sign bit
            cmp AX, 0x401D; // avoid InvalidException in fist
            jae L_extreme;
            fldl2e;
            fmulp ST(1), ST; // y = x*log2(e)
            fist dword ptr [ESP]; // scratchint = rndint(y)
            fisub dword ptr [ESP]; // y - rndint(y)
            // and now set scratchreal exponent
            mov EAX, [ESP];
            add EAX, 0x3fff;
            jle short L_largenegative;
            cmp EAX,0x8000;
            jge short L_largepositive;
            mov [ESP+8+8],AX;
            f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
            fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
            fmul ST(1), ST;  // ST=2rndy, ST(1)=2rndy*2ym1
            fld1;
            fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
            faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
            add ESP,12+8;
            ret PARAMSIZE;

L_extreme:  // Extreme exponent. X is very large positive, very
            // large negative, infinity, or NaN.
            fxam;
            fstsw AX;
            test AX, 0x0400; // NaN_or_zero, but we already know x!=0
            jz L_was_nan;  // if x is NaN, returns x
            test AX, 0x0200;
            jnz L_largenegative;
L_largepositive:
            // Set scratchreal = real.max.
            // squaring it will create infinity, and set overflow flag.
            mov word  ptr [ESP+8+8], 0x7FFE;
            fstp ST(0);
            fld real ptr [ESP+8];  // load scratchreal
            fmul ST(0), ST;        // square it, to create havoc!
L_was_nan:
            add ESP,12+8;
            ret PARAMSIZE;
L_largenegative:
            fstp ST(0);
            fld1;
            fchs; // return -1. Underflow flag is not set.
            add ESP,12+8;
            ret PARAMSIZE;
        }
    }
    else version(D_InlineAsm_X86_64)
    {
        asm
        {
            naked;
        }
        version (Win64)
        {
            asm
            {
                fld   real ptr [RCX];  // x
                mov   AX,[RCX+8];      // AX = exponent and sign
            }
        }
        else
        {
            asm
            {
                fld   real ptr [RSP+8];  // x
                mov   AX,[RSP+8+8];      // AX = exponent and sign
            }
        }
        asm
        {
            /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
             * Author: Don Clugston.
             *
             *    expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
             *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
             *     and 2ym1 = (2^(y-rndint(y))-1).
             *    If 2rndy  < 0.5*real.epsilon, result is -1.
             *    Implementation is otherwise the same as for exp2()
             */
            sub RSP, 24;       // Create scratch space on the stack
            // [RSP,RSP+2] = scratchint
            // [RSP+4..+6, +8..+10, +10] = scratchreal
            // set scratchreal mantissa = 1.0
            mov dword ptr [RSP+8], 0;
            mov dword ptr [RSP+8+4], 0x80000000;
            and AX, 0x7FFF; // drop sign bit
            cmp AX, 0x401D; // avoid InvalidException in fist
            jae L_extreme;
            fldl2e;
            fmul ; // y = x*log2(e)
            fist dword ptr [RSP]; // scratchint = rndint(y)
            fisub dword ptr [RSP]; // y - rndint(y)
            // and now set scratchreal exponent
            mov EAX, [RSP];
            add EAX, 0x3fff;
            jle short L_largenegative;
            cmp EAX,0x8000;
            jge short L_largepositive;
            mov [RSP+8+8],AX;
            f2xm1; // 2^(y-rndint(y)) -1
            fld real ptr [RSP+8] ; // 2^rndint(y)
            fmul ST(1), ST;
            fld1;
            fsubp ST(1), ST;
            fadd;
            add RSP,24;
            ret;

L_extreme:  // Extreme exponent. X is very large positive, very
            // large negative, infinity, or NaN.
            fxam;
            fstsw AX;
            test AX, 0x0400; // NaN_or_zero, but we already know x!=0
            jz L_was_nan;  // if x is NaN, returns x
            test AX, 0x0200;
            jnz L_largenegative;
L_largepositive:
            // Set scratchreal = real.max.
            // squaring it will create infinity, and set overflow flag.
            mov word  ptr [RSP+8+8], 0x7FFE;
            fstp ST(0);
            fld real ptr [RSP+8];  // load scratchreal
            fmul ST(0), ST;        // square it, to create havoc!
L_was_nan:
            add RSP,24;
            ret;

L_largenegative:
            fstp ST(0);
            fld1;
            fchs; // return -1. Underflow flag is not set.
            add RSP,24;
            ret;
        }
    }
    else
    {
        // Coefficients for exp(x) - 1
        static immutable real[5] P = [
           -1.586135578666346600772998894928250240826E4L,
            2.642771505685952966904660652518429479531E3L,
           -3.423199068835684263987132888286791620673E2L,
            1.800826371455042224581246202420972737840E1L,
           -5.238523121205561042771939008061958820811E-1L,
        ];
        static immutable real[6] Q = [
           -9.516813471998079611319047060563358064497E4L,
            3.964866271411091674556850458227710004570E4L,
           -7.207678383830091850230366618190187434796E3L,
            7.206038318724600171970199625081491823079E2L,
           -4.002027679107076077238836622982900945173E1L,
            1.000000000000000000000000000000000000000E0L,
        ];

        // C1 + C2 = LN2.
        enum real C1 = 6.9314575195312500000000E-1L;
        enum real C2 = 1.4286068203094172321215E-6L;

        // Overflow and Underflow limits.
        enum real OF =  1.1356523406294143949492E4L;
        enum real UF = -4.5054566736396445112120088E1L;

        // Special cases.
        if (x > OF)
            return real.infinity;
        if (x == 0.0)
            return x;
        if (x < UF)
            return -1.0;

        // Express x = LN2 (n + remainder), remainder not exceeding 1/2.
        int n = cast(int)floor(0.5 + x / LN2);
        x -= n * C1;
        x -= n * C2;

        // Rational approximation:
        //  exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x)
        real px = x * poly(x, P);
        real qx = poly(x, Q);
        real xx = x * x;
        qx = x + (0.5 * xx + xx * px / qx);

        // We have qx = exp(remainder LN2) - 1, so:
        //  exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1.
        px = ldexp(1.0, n);
        x = px * qx + (px - 1.0);

        return x;
    }
}



/**
 * Calculates 2$(SUP x).
 *
 *  $(TABLE_SV
 *    $(TR $(TH x)             $(TH exp2(x))   )
 *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
 *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
 *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
 *  )
 */
real exp2(real x) @trusted pure nothrow
{
    version(D_InlineAsm_X86)
    {
        enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4

        asm
        {
            /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
             * Author: Don Clugston.
             *
             * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
             * The trick for high performance is to avoid the fscale(28cycles on core2),
             * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
             *
             * We can do frndint by using fist. BUT we can't use it for huge numbers,
             * because it will set the Invalid Operation flag if overflow or NaN occurs.
             * Fortunately, whenever this happens the result would be zero or infinity.
             *
             * We can perform fscale by directly poking into the exponent. BUT this doesn't
             * work for the (very rare) cases where the result is subnormal. So we fall back
             * to the slow method in that case.
             */
            naked;
            fld real ptr [ESP+4] ; // x
            mov AX, [ESP+4+8]; // AX = exponent and sign
            sub ESP, 12+8; // Create scratch space on the stack
            // [ESP,ESP+2] = scratchint
            // [ESP+4..+6, +8..+10, +10] = scratchreal
            // set scratchreal mantissa = 1.0
            mov dword ptr [ESP+8], 0;
            mov dword ptr [ESP+8+4], 0x80000000;
            and AX, 0x7FFF; // drop sign bit
            cmp AX, 0x401D; // avoid InvalidException in fist
            jae L_extreme;
            fist dword ptr [ESP]; // scratchint = rndint(x)
            fisub dword ptr [ESP]; // x - rndint(x)
            // and now set scratchreal exponent
            mov EAX, [ESP];
            add EAX, 0x3fff;
            jle short L_subnormal;
            cmp EAX,0x8000;
            jge short L_overflow;
            mov [ESP+8+8],AX;
L_normal:
            f2xm1;
            fld1;
            faddp ST(1), ST; // 2^^(x-rndint(x))
            fld real ptr [ESP+8] ; // 2^^rndint(x)
            add ESP,12+8;
            fmulp ST(1), ST;
            ret PARAMSIZE;

L_subnormal:
            // Result will be subnormal.
            // In this rare case, the simple poking method doesn't work.
            // The speed doesn't matter, so use the slow fscale method.
            fild dword ptr [ESP];  // scratchint
            fld1;
            fscale;
            fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
            fstp ST(0);         // drop scratchint
            jmp L_normal;

L_extreme:  // Extreme exponent. X is very large positive, very
            // large negative, infinity, or NaN.
            fxam;
            fstsw AX;
            test AX, 0x0400; // NaN_or_zero, but we already know x!=0
            jz L_was_nan;  // if x is NaN, returns x
            // set scratchreal = real.min_normal
            // squaring it will return 0, setting underflow flag
            mov word  ptr [ESP+8+8], 1;
            test AX, 0x0200;
            jnz L_waslargenegative;
L_overflow:
            // Set scratchreal = real.max.
            // squaring it will create infinity, and set overflow flag.
            mov word  ptr [ESP+8+8], 0x7FFE;
L_waslargenegative:
            fstp ST(0);
            fld real ptr [ESP+8];  // load scratchreal
            fmul ST(0), ST;        // square it, to create havoc!
L_was_nan:
            add ESP,12+8;
            ret PARAMSIZE;
        }
    }
    else version(D_InlineAsm_X86_64)
    {
        asm
        {
            naked;
        }
        version (Win64)
        {
            asm
            {
                fld   real ptr [RCX];  // x
                mov   AX,[RCX+8];      // AX = exponent and sign
            }
        }
        else
        {
            asm
            {
                fld   real ptr [RSP+8];  // x
                mov   AX,[RSP+8+8];      // AX = exponent and sign
            }
        }
        asm
        {
            /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
             * Author: Don Clugston.
             *
             * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
             * The trick for high performance is to avoid the fscale(28cycles on core2),
             * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
             *
             * We can do frndint by using fist. BUT we can't use it for huge numbers,
             * because it will set the Invalid Operation flag is overflow or NaN occurs.
             * Fortunately, whenever this happens the result would be zero or infinity.
             *
             * We can perform fscale by directly poking into the exponent. BUT this doesn't
             * work f…
Summary ✨

This D code implements mathematical functions and checks for numerical precision, including yl2x and yl2xp1, which translate to a single x87 instruction, and various tests for floating-point operations such as sqrt, fabs, sin, cos, and tan. It also includes tests for the floor and ceil functions in pure functions.
Alerts (6)

Complexity hotspot; lines 334 to 339 (total complexity: 6)
334 335 336 337 338 339