PageRenderTime 147ms CodeModel.GetById 4ms app.highlight 124ms RepoModel.GetById 2ms app.codeStats 1ms

/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 files are truncated, but you can click here to view the full file

   1// Written in the D programming language.
   2
   3/**
   4 * Elementary mathematical functions
   5 *
   6 * Contains the elementary mathematical functions (powers, roots,
   7 * and trignometric functions), and low-level floating-point operations.
   8 * Mathematical special functions are available in std.mathspecial.
   9 *
  10 * The functionality closely follows the IEEE754-2008 standard for
  11 * floating-point arithmetic, including the use of camelCase names rather
  12 * than C99-style lower case names. All of these functions behave correctly
  13 * when presented with an infinity or NaN.
  14 *
  15 * Unlike C, there is no global 'errno' variable. Consequently, almost all of
  16 * these functions are pure nothrow.
  17 *
  18 * Status:
  19 * The semantics and names of feqrel and approxEqual will be revised.
  20 *
  21 * Macros:
  22 *      WIKI = Phobos/StdMath
  23 *
  24 *      TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
  25 *              <caption>Special Values</caption>
  26 *              $0</table>
  27 *      SVH = $(TR $(TH $1) $(TH $2))
  28 *      SV  = $(TR $(TD $1) $(TD $2))
  29 *
  30 *      NAN = $(RED NAN)
  31 *      SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
  32 *      GAMMA = &#915;
  33 *      THETA = &theta;
  34 *      INTEGRAL = &#8747;
  35 *      INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
  36 *      POWER = $1<sup>$2</sup>
  37 *      SUB = $1<sub>$2</sub>
  38 *      BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
  39 *      CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
  40 *      PLUSMN = &plusmn;
  41 *      INFIN = &infin;
  42 *      PLUSMNINF = &plusmn;&infin;
  43 *      PI = &pi;
  44 *      LT = &lt;
  45 *      GT = &gt;
  46 *      SQRT = &radic;
  47 *      HALF = &frac12;
  48 *
  49 * Copyright: Copyright Digital Mars 2000 - 2011.
  50 *            D implementations of tan, atan, atan2, exp, expm1, exp2, log, log10, log1p,
  51 *            log2, floor, ceil and lrint functions are based on the CEPHES math library,
  52 *            which is Copyright (C) 2001 Stephen L. Moshier <steve@moshier.net>
  53 *            and are incorporated herein by permission of the author.  The author
  54 *            reserves the right to distribute this material elsewhere under different
  55 *            copying permissions.  These modifications are distributed here under
  56 *            the following terms:
  57 * License:   <a href="http://www.boost.org/LICENSE_1_0.txt">Boost License 1.0</a>.
  58 * Authors:   $(WEB digitalmars.com, Walter Bright),
  59 *                        Don Clugston, Conversion of CEPHES math library to D by Iain Buclaw
  60 * Source: $(PHOBOSSRC std/_math.d)
  61 */
  62module std.math;
  63
  64import core.stdc.math;
  65import std.traits;
  66
  67version(unittest)
  68{
  69    import std.typetuple;
  70}
  71
  72version(LDC)
  73{
  74    import ldc.intrinsics;
  75}
  76
  77version(DigitalMars)
  78{
  79    version = INLINE_YL2X;        // x87 has opcodes for these
  80}
  81
  82version (X86)
  83{
  84    version = X86_Any;
  85}
  86
  87version (X86_64)
  88{
  89    version = X86_Any;
  90}
  91
  92version(D_InlineAsm_X86)
  93{
  94    version = InlineAsm_X86_Any;
  95}
  96else version(D_InlineAsm_X86_64)
  97{
  98    version = InlineAsm_X86_Any;
  99}
 100
 101
 102version(unittest)
 103{
 104    import core.stdc.stdio;
 105
 106    static if(real.sizeof > double.sizeof)
 107        enum uint useDigits = 16;
 108    else
 109        enum uint useDigits = 15;
 110
 111    /******************************************
 112     * Compare floating point numbers to n decimal digits of precision.
 113     * Returns:
 114     *  1       match
 115     *  0       nomatch
 116     */
 117
 118    private bool equalsDigit(real x, real y, uint ndigits)
 119    {
 120        if (signbit(x) != signbit(y))
 121            return 0;
 122
 123        if (isinf(x) && isinf(y))
 124            return 1;
 125        if (isinf(x) || isinf(y))
 126            return 0;
 127
 128        if (isnan(x) && isnan(y))
 129            return 1;
 130        if (isnan(x) || isnan(y))
 131            return 0;
 132
 133        char[30] bufx;
 134        char[30] bufy;
 135        assert(ndigits < bufx.length);
 136
 137        int ix;
 138        int iy;
 139        ix = sprintf(bufx.ptr, "%.*Lg", ndigits, x);
 140        assert(ix < bufx.length && ix > 0);
 141        iy = sprintf(bufy.ptr, "%.*Lg", ndigits, y);
 142        assert(ix < bufy.length && ix > 0);
 143
 144        return bufx[0 .. ix] == bufy[0 .. iy];
 145    }
 146}
 147
 148
 149
 150private:
 151/*
 152 * The following IEEE 'real' formats are currently supported:
 153 * 64 bit Big-endian  'double' (eg PowerPC)
 154 * 128 bit Big-endian 'quadruple' (eg SPARC)
 155 * 64 bit Little-endian 'double' (eg x86-SSE2)
 156 * 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium).
 157 * 128 bit Little-endian 'quadruple' (not implemented on any known processor!)
 158 *
 159 * Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support
 160 */
 161version(LittleEndian)
 162{
 163    static assert(real.mant_dig == 53 || real.mant_dig==64
 164               || real.mant_dig == 113,
 165      "Only 64-bit, 80-bit, and 128-bit reals"
 166      " are supported for LittleEndian CPUs");
 167}
 168else
 169{
 170    static assert(real.mant_dig == 53 || real.mant_dig==106
 171               || real.mant_dig == 113,
 172    "Only 64-bit and 128-bit reals are supported for BigEndian CPUs."
 173    " double-double reals have partial support");
 174}
 175
 176// Constants used for extracting the components of the representation.
 177// They supplement the built-in floating point properties.
 178template floatTraits(T)
 179{
 180    // EXPMASK is a ushort mask to select the exponent portion (without sign)
 181    // EXPPOS_SHORT is the index of the exponent when represented as a ushort array.
 182    // SIGNPOS_BYTE is the index of the sign when represented as a ubyte array.
 183    // RECIP_EPSILON is the value such that (smallest_subnormal) * RECIP_EPSILON == T.min_normal
 184    enum T RECIP_EPSILON = (1/T.epsilon);
 185    static if (T.mant_dig == 24)
 186    { // float
 187        enum ushort EXPMASK = 0x7F80;
 188        enum ushort EXPBIAS = 0x3F00;
 189        enum uint EXPMASK_INT = 0x7F80_0000;
 190        enum uint MANTISSAMASK_INT = 0x007F_FFFF;
 191        version(LittleEndian)
 192        {
 193            enum EXPPOS_SHORT = 1;
 194        }
 195        else
 196        {
 197            enum EXPPOS_SHORT = 0;
 198        }
 199    }
 200    else static if (T.mant_dig == 53) // double, or real==double
 201    {
 202        enum ushort EXPMASK = 0x7FF0;
 203        enum ushort EXPBIAS = 0x3FE0;
 204        enum uint EXPMASK_INT = 0x7FF0_0000;
 205        enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only
 206        version(LittleEndian)
 207        {
 208            enum EXPPOS_SHORT = 3;
 209            enum SIGNPOS_BYTE = 7;
 210        }
 211        else
 212        {
 213            enum EXPPOS_SHORT = 0;
 214            enum SIGNPOS_BYTE = 0;
 215        }
 216    }
 217    else static if (T.mant_dig == 64) // real80
 218    {
 219        enum ushort EXPMASK = 0x7FFF;
 220        enum ushort EXPBIAS = 0x3FFE;
 221        version(LittleEndian)
 222        {
 223            enum EXPPOS_SHORT = 4;
 224            enum SIGNPOS_BYTE = 9;
 225        }
 226        else
 227        {
 228            enum EXPPOS_SHORT = 0;
 229            enum SIGNPOS_BYTE = 0;
 230        }
 231    }
 232    else static if (T.mant_dig == 113) // quadruple
 233    {
 234        enum ushort EXPMASK = 0x7FFF;
 235        version(LittleEndian)
 236        {
 237            enum EXPPOS_SHORT = 7;
 238            enum SIGNPOS_BYTE = 15;
 239        }
 240        else
 241        {
 242            enum EXPPOS_SHORT = 0;
 243            enum SIGNPOS_BYTE = 0;
 244        }
 245    }
 246    else static if (T.mant_dig == 106) // doubledouble
 247    {
 248        enum ushort EXPMASK = 0x7FF0;
 249        // the exponent byte is not unique
 250        version(LittleEndian)
 251        {
 252            enum EXPPOS_SHORT = 7; // [3] is also an exp short
 253            enum SIGNPOS_BYTE = 15;
 254        }
 255        else
 256        {
 257            enum EXPPOS_SHORT = 0; // [4] is also an exp short
 258            enum SIGNPOS_BYTE = 0;
 259        }
 260    }
 261}
 262
 263// These apply to all floating-point types
 264version(LittleEndian)
 265{
 266    enum MANTISSA_LSB = 0;
 267    enum MANTISSA_MSB = 1;
 268}
 269else
 270{
 271    enum MANTISSA_LSB = 1;
 272    enum MANTISSA_MSB = 0;
 273}
 274
 275public:
 276
 277// Values obtained from Wolfram Alpha. 116 bits ought to be enough for anybody.
 278// Wolfram Alpha LLC. 2011. Wolfram|Alpha. http://www.wolframalpha.com/input/?i=e+in+base+16 (access July 6, 2011).
 279enum real E =          0x1.5bf0a8b1457695355fb8ac404e7a8p+1L; /** e = 2.718281... */
 280enum real LOG2T =      0x1.a934f0979a3715fc9257edfe9b5fbp+1L; /** $(SUB log, 2)10 = 3.321928... */
 281enum real LOG2E =      0x1.71547652b82fe1777d0ffda0d23a8p+0L; /** $(SUB log, 2)e = 1.442695... */
 282enum real LOG2 =       0x1.34413509f79fef311f12b35816f92p-2L; /** $(SUB log, 10)2 = 0.301029... */
 283enum real LOG10E =     0x1.bcb7b1526e50e32a6ab7555f5a67cp-2L; /** $(SUB log, 10)e = 0.434294... */
 284enum real LN2 =        0x1.62e42fefa39ef35793c7673007e5fp-1L; /** ln 2  = 0.693147... */
 285enum real LN10 =       0x1.26bb1bbb5551582dd4adac5705a61p+1L; /** ln 10 = 2.302585... */
 286enum real PI =         0x1.921fb54442d18469898cc51701b84p+1L; /** $(_PI) = 3.141592... */
 287enum real PI_2 =       PI/2;                                  /** $(PI) / 2 = 1.570796... */
 288enum real PI_4 =       PI/4;                                  /** $(PI) / 4 = 0.785398... */
 289enum real M_1_PI =     0x1.45f306dc9c882a53f84eafa3ea69cp-2L; /** 1 / $(PI) = 0.318309... */
 290enum real M_2_PI =     2*M_1_PI;                              /** 2 / $(PI) = 0.636619... */
 291enum real M_2_SQRTPI = 0x1.20dd750429b6d11ae3a914fed7fd8p+0L; /** 2 / $(SQRT)$(PI) = 1.128379... */
 292enum real SQRT2 =      0x1.6a09e667f3bcc908b2fb1366ea958p+0L; /** $(SQRT)2 = 1.414213... */
 293enum real SQRT1_2 =    SQRT2/2;                               /** $(SQRT)$(HALF) = 0.707106... */
 294// Note: Make sure the magic numbers in compiler backend for x87 match these.
 295
 296
 297/***********************************
 298 * Calculates the absolute value
 299 *
 300 * For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) )
 301 * = hypot(z.re, z.im).
 302 */
 303Num abs(Num)(Num x) @safe pure nothrow
 304    if (is(typeof(Num.init >= 0)) && is(typeof(-Num.init)) &&
 305            !(is(Num* : const(ifloat*)) || is(Num* : const(idouble*))
 306                    || is(Num* : const(ireal*))))
 307{
 308    static if (isFloatingPoint!(Num))
 309        return fabs(x);
 310    else
 311        return x>=0 ? x : -x;
 312}
 313
 314auto abs(Num)(Num z) @safe pure nothrow
 315    if (is(Num* : const(cfloat*)) || is(Num* : const(cdouble*))
 316            || is(Num* : const(creal*)))
 317{
 318    return hypot(z.re, z.im);
 319}
 320
 321/** ditto */
 322real abs(Num)(Num y) @safe pure nothrow
 323    if (is(Num* : const(ifloat*)) || is(Num* : const(idouble*))
 324            || is(Num* : const(ireal*)))
 325{
 326    return fabs(y.im);
 327}
 328
 329
 330unittest
 331{
 332    assert(isIdentical(abs(-0.0L), 0.0L));
 333    assert(isNaN(abs(real.nan)));
 334    assert(abs(-real.infinity) == real.infinity);
 335    assert(abs(-3.2Li) == 3.2L);
 336    assert(abs(71.6Li) == 71.6L);
 337    assert(abs(-56) == 56);
 338    assert(abs(2321312L)  == 2321312L);
 339    assert(abs(-1+1i) == sqrt(2.0L));
 340}
 341
 342/***********************************
 343 * Complex conjugate
 344 *
 345 *  conj(x + iy) = x - iy
 346 *
 347 * Note that z * conj(z) = $(POWER z.re, 2) - $(POWER z.im, 2)
 348 * is always a real number
 349 */
 350creal conj(creal z) @safe pure nothrow
 351{
 352    return z.re - z.im*1i;
 353}
 354
 355/** ditto */
 356ireal conj(ireal y) @safe pure nothrow
 357{
 358    return -y;
 359}
 360
 361
 362unittest
 363{
 364    assert(conj(7 + 3i) == 7-3i);
 365    ireal z = -3.2Li;
 366    assert(conj(z) == -z);
 367}
 368
 369/***********************************
 370 * Returns cosine of x. x is in radians.
 371 *
 372 *      $(TABLE_SV
 373 *      $(TR $(TH x)                 $(TH cos(x)) $(TH invalid?))
 374 *      $(TR $(TD $(NAN))            $(TD $(NAN)) $(TD yes)     )
 375 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)     )
 376 *      )
 377 * Bugs:
 378 *      Results are undefined if |x| >= $(POWER 2,64).
 379 */
 380
 381real cos(real x) @safe pure nothrow;       /* intrinsic */
 382
 383/***********************************
 384 * Returns sine of x. x is in radians.
 385 *
 386 *      $(TABLE_SV
 387 *      $(TR $(TH x)               $(TH sin(x))      $(TH invalid?))
 388 *      $(TR $(TD $(NAN))          $(TD $(NAN))      $(TD yes))
 389 *      $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0) $(TD no))
 390 *      $(TR $(TD $(PLUSMNINF))    $(TD $(NAN))      $(TD yes))
 391 *      )
 392 * Bugs:
 393 *      Results are undefined if |x| >= $(POWER 2,64).
 394 */
 395
 396real sin(real x) @safe pure nothrow;       /* intrinsic */
 397
 398
 399/***********************************
 400 *  sine, complex and imaginary
 401 *
 402 *  sin(z) = sin(z.re)*cosh(z.im) + cos(z.re)*sinh(z.im)i
 403 *
 404 * If both sin($(THETA)) and cos($(THETA)) are required,
 405 * it is most efficient to use expi($(THETA)).
 406 */
 407creal sin(creal z) @safe pure nothrow
 408{
 409    creal cs = expi(z.re);
 410    creal csh = coshisinh(z.im);
 411    return cs.im * csh.re + cs.re * csh.im * 1i;
 412}
 413
 414/** ditto */
 415ireal sin(ireal y) @safe pure nothrow
 416{
 417    return cosh(y.im)*1i;
 418}
 419
 420unittest
 421{
 422  assert(sin(0.0+0.0i) == 0.0);
 423  assert(sin(2.0+0.0i) == sin(2.0L) );
 424}
 425
 426/***********************************
 427 *  cosine, complex and imaginary
 428 *
 429 *  cos(z) = cos(z.re)*cosh(z.im) - sin(z.re)*sinh(z.im)i
 430 */
 431creal cos(creal z) @safe pure nothrow
 432{
 433    creal cs = expi(z.re);
 434    creal csh = coshisinh(z.im);
 435    return cs.re * csh.re - cs.im * csh.im * 1i;
 436}
 437
 438/** ditto */
 439real cos(ireal y) @safe pure nothrow
 440{
 441    return cosh(y.im);
 442}
 443
 444unittest
 445{
 446    assert(cos(0.0+0.0i)==1.0);
 447    assert(cos(1.3L+0.0i)==cos(1.3L));
 448    assert(cos(5.2Li)== cosh(5.2L));
 449}
 450
 451/****************************************************************************
 452 * Returns tangent of x. x is in radians.
 453 *
 454 *      $(TABLE_SV
 455 *      $(TR $(TH x)             $(TH tan(x))       $(TH invalid?))
 456 *      $(TR $(TD $(NAN))        $(TD $(NAN))       $(TD yes))
 457 *      $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) $(TD no))
 458 *      $(TR $(TD $(PLUSMNINF))  $(TD $(NAN))       $(TD yes))
 459 *      )
 460 */
 461
 462real tan(real x) @trusted pure nothrow
 463{
 464    version(D_InlineAsm_X86)
 465    {
 466    asm
 467    {
 468        fld     x[EBP]                  ; // load theta
 469        fxam                            ; // test for oddball values
 470        fstsw   AX                      ;
 471        sahf                            ;
 472        jc      trigerr                 ; // x is NAN, infinity, or empty
 473                                          // 387's can handle subnormals
 474SC18:   fptan                           ;
 475        fstp    ST(0)                   ; // dump X, which is always 1
 476        fstsw   AX                      ;
 477        sahf                            ;
 478        jnp     Lret                    ; // C2 = 1 (x is out of range)
 479
 480        // Do argument reduction to bring x into range
 481        fldpi                           ;
 482        fxch                            ;
 483SC17:   fprem1                          ;
 484        fstsw   AX                      ;
 485        sahf                            ;
 486        jp      SC17                    ;
 487        fstp    ST(1)                   ; // remove pi from stack
 488        jmp     SC18                    ;
 489
 490trigerr:
 491        jnp     Lret                    ; // if theta is NAN, return theta
 492        fstp    ST(0)                   ; // dump theta
 493    }
 494    return real.nan;
 495
 496Lret: {}
 497    }
 498    else version(D_InlineAsm_X86_64)
 499    {
 500        version (Win64)
 501        {
 502            asm
 503            {
 504                fld     real ptr [RCX]  ; // load theta
 505            }
 506        }
 507        else
 508        {
 509            asm
 510            {
 511                fld     x[RBP]          ; // load theta
 512            }
 513        }
 514    asm
 515    {
 516        fxam                            ; // test for oddball values
 517        fstsw   AX                      ;
 518        test    AH,1                    ;
 519        jnz     trigerr                 ; // x is NAN, infinity, or empty
 520                                          // 387's can handle subnormals
 521SC18:   fptan                           ;
 522        fstp    ST(0)                   ; // dump X, which is always 1
 523        fstsw   AX                      ;
 524        test    AH,4                    ;
 525        jz      Lret                    ; // C2 = 1 (x is out of range)
 526
 527        // Do argument reduction to bring x into range
 528        fldpi                           ;
 529        fxch                            ;
 530SC17:   fprem1                          ;
 531        fstsw   AX                      ;
 532        test    AH,4                    ;
 533        jnz     SC17                    ;
 534        fstp    ST(1)                   ; // remove pi from stack
 535        jmp     SC18                    ;
 536
 537trigerr:
 538        test    AH,4                    ;
 539        jz      Lret                    ; // if theta is NAN, return theta
 540        fstp    ST(0)                   ; // dump theta
 541    }
 542    return real.nan;
 543
 544Lret: {}
 545    }
 546    else
 547    {
 548        // Coefficients for tan(x)
 549        static immutable real[3] P = [
 550           -1.7956525197648487798769E7L,
 551            1.1535166483858741613983E6L,
 552           -1.3093693918138377764608E4L,
 553        ];
 554        static immutable real[5] Q = [
 555           -5.3869575592945462988123E7L,
 556            2.5008380182335791583922E7L,
 557           -1.3208923444021096744731E6L,
 558            1.3681296347069295467845E4L,
 559            1.0000000000000000000000E0L,
 560        ];
 561
 562        // PI/4 split into three parts.
 563        enum real P1 = 7.853981554508209228515625E-1L;
 564        enum real P2 = 7.946627356147928367136046290398E-9L;
 565        enum real P3 = 3.061616997868382943065164830688E-17L;
 566
 567        // Special cases.
 568        if (x == 0.0 || isNaN(x))
 569            return x;
 570        if (isInfinity(x))
 571            return real.nan;
 572
 573        // Make argument positive but save the sign.
 574        bool sign = false;
 575        if (signbit(x))
 576        {
 577            sign = true;
 578            x = -x;
 579        }
 580
 581        // Compute x mod PI/4.
 582        real y = floor(x / PI_4);
 583        // Strip high bits of integer part.
 584        real z = ldexp(y, -4);
 585        // Compute y - 16 * (y / 16).
 586        z = y - ldexp(floor(z), 4);
 587
 588        // Integer and fraction part modulo one octant.
 589        int j = cast(int)(z);
 590
 591        // Map zeros and singularities to origin.
 592        if (j & 1)
 593        {
 594            j += 1;
 595            y += 1.0;
 596        }
 597
 598        z = ((x - y * P1) - y * P2) - y * P3;
 599        real zz = z * z;
 600
 601        if (zz > 1.0e-20L)
 602            y = z + z * (zz * poly(zz, P) / poly(zz, Q));
 603        else
 604            y = z;
 605
 606        if (j & 2)
 607            y = -1.0 / y;
 608
 609        return (sign) ? -y : y;
 610    }
 611}
 612
 613unittest
 614{
 615    static real[2][] vals =     // angle,tan
 616        [
 617         [   0,   0],
 618         [   .5,  .5463024898],
 619         [   1,   1.557407725],
 620         [   1.5, 14.10141995],
 621         [   2,  -2.185039863],
 622         [   2.5,-.7470222972],
 623         [   3,  -.1425465431],
 624         [   3.5, .3745856402],
 625         [   4,   1.157821282],
 626         [   4.5, 4.637332055],
 627         [   5,  -3.380515006],
 628         [   5.5,-.9955840522],
 629         [   6,  -.2910061914],
 630         [   6.5, .2202772003],
 631         [   10,  .6483608275],
 632
 633         // special angles
 634         [   PI_4,   1],
 635         //[   PI_2,   real.infinity], // PI_2 is not _exactly_ pi/2.
 636         [   3*PI_4, -1],
 637         [   PI,     0],
 638         [   5*PI_4, 1],
 639         //[   3*PI_2, -real.infinity],
 640         [   7*PI_4, -1],
 641         [   2*PI,   0],
 642         ];
 643    int i;
 644
 645    for (i = 0; i < vals.length; i++)
 646    {
 647        real x = vals[i][0];
 648        real r = vals[i][1];
 649        real t = tan(x);
 650
 651        //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
 652        if (!isIdentical(r, t)) assert(fabs(r-t) <= .0000001);
 653
 654        x = -x;
 655        r = -r;
 656        t = tan(x);
 657        //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
 658        if (!isIdentical(r, t) && !(r!=r && t!=t)) assert(fabs(r-t) <= .0000001);
 659    }
 660    // overflow
 661    assert(isNaN(tan(real.infinity)));
 662    assert(isNaN(tan(-real.infinity)));
 663    // NaN propagation
 664    assert(isIdentical( tan(NaN(0x0123L)), NaN(0x0123L) ));
 665}
 666
 667unittest
 668{
 669    assert(equalsDigit(tan(PI / 3), std.math.sqrt(3.0), useDigits));
 670}
 671
 672/***************
 673 * Calculates the arc cosine of x,
 674 * returning a value ranging from 0 to $(PI).
 675 *
 676 *      $(TABLE_SV
 677 *      $(TR $(TH x)         $(TH acos(x)) $(TH invalid?))
 678 *      $(TR $(TD $(GT)1.0)  $(TD $(NAN))  $(TD yes))
 679 *      $(TR $(TD $(LT)-1.0) $(TD $(NAN))  $(TD yes))
 680 *      $(TR $(TD $(NAN))    $(TD $(NAN))  $(TD yes))
 681 *  )
 682 */
 683real acos(real x) @safe pure nothrow
 684{
 685    return atan2(sqrt(1-x*x), x);
 686}
 687
 688/// ditto
 689double acos(double x) @safe pure nothrow { return acos(cast(real)x); }
 690
 691/// ditto
 692float acos(float x) @safe pure nothrow  { return acos(cast(real)x); }
 693
 694unittest
 695{
 696    assert(equalsDigit(acos(0.5), std.math.PI / 3, useDigits));
 697}
 698
 699/***************
 700 * Calculates the arc sine of x,
 701 * returning a value ranging from -$(PI)/2 to $(PI)/2.
 702 *
 703 *      $(TABLE_SV
 704 *      $(TR $(TH x)            $(TH asin(x))      $(TH invalid?))
 705 *      $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
 706 *      $(TR $(TD $(GT)1.0)     $(TD $(NAN))       $(TD yes))
 707 *      $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD yes))
 708 *  )
 709 */
 710real asin(real x) @safe pure nothrow
 711{
 712    return atan2(x, sqrt(1-x*x));
 713}
 714
 715/// ditto
 716double asin(double x) @safe pure nothrow { return asin(cast(real)x); }
 717
 718/// ditto
 719float asin(float x) @safe pure nothrow  { return asin(cast(real)x); }
 720
 721unittest
 722{
 723    assert(equalsDigit(asin(0.5), PI / 6, useDigits));
 724}
 725
 726/***************
 727 * Calculates the arc tangent of x,
 728 * returning a value ranging from -$(PI)/2 to $(PI)/2.
 729 *
 730 *      $(TABLE_SV
 731 *      $(TR $(TH x)                 $(TH atan(x))      $(TH invalid?))
 732 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no))
 733 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))       $(TD yes))
 734 *  )
 735 */
 736real atan(real x) @safe pure nothrow
 737{
 738    version(InlineAsm_X86_Any)
 739    {
 740        return atan2(x, 1.0L);
 741    }
 742    else
 743    {
 744        // Coefficients for atan(x)
 745        static immutable real[5] P = [
 746           -5.0894116899623603312185E1L,
 747           -9.9988763777265819915721E1L,
 748           -6.3976888655834347413154E1L,
 749           -1.4683508633175792446076E1L,
 750           -8.6863818178092187535440E-1L,
 751        ];
 752        static immutable real[6] Q = [
 753            1.5268235069887081006606E2L,
 754            3.9157570175111990631099E2L,
 755            3.6144079386152023162701E2L,
 756            1.4399096122250781605352E2L,
 757            2.2981886733594175366172E1L,
 758            1.0000000000000000000000E0L,
 759        ];
 760
 761        // tan(PI/8)
 762        enum real TAN_PI_8 = 4.1421356237309504880169e-1L;
 763        // tan(3 * PI/8)
 764        enum real TAN3_PI_8 = 2.41421356237309504880169L;
 765
 766        // Special cases.
 767        if (x == 0.0)
 768            return x;
 769        if (isInfinity(x))
 770            return copysign(PI_2, x);
 771
 772        // Make argument positive but save the sign.
 773        bool sign = false;
 774        if (signbit(x))
 775        {
 776            sign = true;
 777            x = -x;
 778        }
 779
 780        // Range reduction.
 781        real y;
 782        if (x > TAN3_PI_8)
 783        {
 784            y = PI_2;
 785            x = -(1.0 / x);
 786        }
 787        else if (x > TAN_PI_8)
 788        {
 789            y = PI_4;
 790            x = (x - 1.0)/(x + 1.0);
 791        }
 792        else
 793            y = 0.0;
 794
 795        // Rational form in x^^2.
 796        real z = x * x;
 797        y = y + (poly(z, P) / poly(z, Q)) * z * x + x;
 798
 799        return (sign) ? -y : y;
 800    }
 801}
 802
 803/// ditto
 804double atan(double x) @safe pure nothrow { return atan(cast(real)x); }
 805
 806/// ditto
 807float atan(float x)  @safe pure nothrow { return atan(cast(real)x); }
 808
 809unittest
 810{
 811    assert(equalsDigit(atan(std.math.sqrt(3.0)), PI / 3, useDigits));
 812}
 813
 814/***************
 815 * Calculates the arc tangent of y / x,
 816 * returning a value ranging from -$(PI) to $(PI).
 817 *
 818 *      $(TABLE_SV
 819 *      $(TR $(TH y)                 $(TH x)            $(TH atan(y, x)))
 820 *      $(TR $(TD $(NAN))            $(TD anything)     $(TD $(NAN)) )
 821 *      $(TR $(TD anything)          $(TD $(NAN))       $(TD $(NAN)) )
 822 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(GT)0.0)     $(TD $(PLUSMN)0.0) )
 823 *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0)         $(TD $(PLUSMN)0.0) )
 824 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(LT)0.0)     $(TD $(PLUSMN)$(PI)))
 825 *      $(TR $(TD $(PLUSMN)0.0)      $(TD -0.0)         $(TD $(PLUSMN)$(PI)))
 826 *      $(TR $(TD $(GT)0.0)          $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
 827 *      $(TR $(TD $(LT)0.0)          $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) )
 828 *      $(TR $(TD $(GT)0.0)          $(TD $(INFIN))     $(TD $(PLUSMN)0.0) )
 829 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything)     $(TD $(PLUSMN)$(PI)/2))
 830 *      $(TR $(TD $(GT)0.0)          $(TD -$(INFIN))    $(TD $(PLUSMN)$(PI)) )
 831 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN))     $(TD $(PLUSMN)$(PI)/4))
 832 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN))    $(TD $(PLUSMN)3$(PI)/4))
 833 *      )
 834 */
 835real atan2(real y, real x) @trusted pure nothrow
 836{
 837    version(InlineAsm_X86_Any)
 838    {
 839        version (Win64)
 840        {
 841            asm {
 842                naked;
 843                fld real ptr [RDX]; // y
 844                fld real ptr [RCX]; // x
 845                fpatan;
 846                ret;
 847            }
 848        }
 849        else
 850        {
 851            asm {
 852                fld y;
 853                fld x;
 854                fpatan;
 855            }
 856        }
 857    }
 858    else
 859    {
 860        // Special cases.
 861        if (isNaN(x) || isNaN(y))
 862            return real.nan;
 863        if (y == 0.0)
 864        {
 865            if (x >= 0 && !signbit(x))
 866                return copysign(0, y);
 867            else
 868                return copysign(PI, y);
 869        }
 870        if (x == 0.0)
 871            return copysign(PI_2, y);
 872        if (isInfinity(x))
 873        {
 874            if (signbit(x))
 875            {
 876                if (isInfinity(y))
 877                    return copysign(3*PI_4, y);
 878                else
 879                    return copysign(PI, y);
 880            }
 881            else
 882            {
 883                if (isInfinity(y))
 884                    return copysign(PI_4, y);
 885                else
 886                    return copysign(0.0, y);
 887            }
 888        }
 889        if (isInfinity(y))
 890            return copysign(PI_2, y);
 891
 892        // Call atan and determine the quadrant.
 893        real z = atan(y / x);
 894
 895        if (signbit(x))
 896        {
 897            if (signbit(y))
 898                z = z - PI;
 899            else
 900                z = z + PI;
 901        }
 902
 903        if (z == 0.0)
 904            return copysign(z, y);
 905
 906        return z;
 907    }
 908}
 909
 910/// ditto
 911double atan2(double y, double x) @safe pure nothrow
 912{
 913    return atan2(cast(real)y, cast(real)x);
 914}
 915
 916/// ditto
 917float atan2(float y, float x) @safe pure nothrow
 918{
 919    return atan2(cast(real)y, cast(real)x);
 920}
 921
 922unittest
 923{
 924    assert(equalsDigit(atan2(1.0L, std.math.sqrt(3.0L)), PI / 6, useDigits));
 925}
 926
 927/***********************************
 928 * Calculates the hyperbolic cosine of x.
 929 *
 930 *      $(TABLE_SV
 931 *      $(TR $(TH x)                 $(TH cosh(x))      $(TH invalid?))
 932 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
 933 *      )
 934 */
 935real cosh(real x) @safe pure nothrow
 936{
 937    //  cosh = (exp(x)+exp(-x))/2.
 938    // The naive implementation works correctly.
 939    real y = exp(x);
 940    return (y + 1.0/y) * 0.5;
 941}
 942
 943/// ditto
 944double cosh(double x) @safe pure nothrow { return cosh(cast(real)x); }
 945
 946/// ditto
 947float cosh(float x) @safe pure nothrow  { return cosh(cast(real)x); }
 948
 949unittest
 950{
 951    assert(equalsDigit(cosh(1.0), (E + 1.0 / E) / 2, useDigits));
 952}
 953
 954/***********************************
 955 * Calculates the hyperbolic sine of x.
 956 *
 957 *      $(TABLE_SV
 958 *      $(TR $(TH x)                 $(TH sinh(x))           $(TH invalid?))
 959 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no))
 960 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
 961 *      )
 962 */
 963real sinh(real x) @safe pure nothrow
 964{
 965    //  sinh(x) =  (exp(x)-exp(-x))/2;
 966    // Very large arguments could cause an overflow, but
 967    // the maximum value of x for which exp(x) + exp(-x)) != exp(x)
 968    // is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80.
 969    if (fabs(x) > real.mant_dig * LN2)
 970    {
 971        return copysign(0.5 * exp(fabs(x)), x);
 972    }
 973
 974    real y = expm1(x);
 975    return 0.5 * y / (y+1) * (y+2);
 976}
 977
 978/// ditto
 979double sinh(double x) @safe pure nothrow { return sinh(cast(real)x); }
 980
 981/// ditto
 982float sinh(float x) @safe pure nothrow  { return sinh(cast(real)x); }
 983
 984unittest
 985{
 986    assert(equalsDigit(sinh(1.0), (E - 1.0 / E) / 2, useDigits));
 987}
 988
 989/***********************************
 990 * Calculates the hyperbolic tangent of x.
 991 *
 992 *      $(TABLE_SV
 993 *      $(TR $(TH x)                 $(TH tanh(x))      $(TH invalid?))
 994 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no) )
 995 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
 996 *      )
 997 */
 998real tanh(real x) @safe pure nothrow
 999{
1000    //  tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
1001    if (fabs(x) > real.mant_dig * LN2)
1002    {
1003        return copysign(1, x);
1004    }
1005
1006    real y = expm1(2*x);
1007    return y / (y + 2);
1008}
1009
1010/// ditto
1011double tanh(double x) @safe pure nothrow { return tanh(cast(real)x); }
1012
1013/// ditto
1014float tanh(float x) @safe pure nothrow { return tanh(cast(real)x); }
1015
1016unittest
1017{
1018    assert(equalsDigit(tanh(1.0), sinh(1.0) / cosh(1.0), 15));
1019}
1020
1021package:
1022
1023/* Returns cosh(x) + I * sinh(x)
1024 * Only one call to exp() is performed.
1025 */
1026creal coshisinh(real x) @safe pure nothrow
1027{
1028    // See comments for cosh, sinh.
1029    if (fabs(x) > real.mant_dig * LN2)
1030    {
1031        real y = exp(fabs(x));
1032        return y * 0.5 + 0.5i * copysign(y, x);
1033    }
1034    else
1035    {
1036        real y = expm1(x);
1037        return (y + 1.0 + 1.0/(y + 1.0)) * 0.5 + 0.5i * y / (y+1) * (y+2);
1038    }
1039}
1040
1041unittest
1042{
1043    creal c = coshisinh(3.0L);
1044    assert(c.re == cosh(3.0L));
1045    assert(c.im == sinh(3.0L));
1046}
1047
1048public:
1049
1050/***********************************
1051 * Calculates the inverse hyperbolic cosine of x.
1052 *
1053 *  Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
1054 *
1055 * $(TABLE_DOMRG
1056 *  $(DOMAIN 1..$(INFIN))
1057 *  $(RANGE  1..log(real.max), $(INFIN)) )
1058 *      $(TABLE_SV
1059 *    $(SVH  x,     acosh(x) )
1060 *    $(SV  $(NAN), $(NAN) )
1061 *    $(SV  $(LT)1,     $(NAN) )
1062 *    $(SV  1,      0       )
1063 *    $(SV  +$(INFIN),+$(INFIN))
1064 *  )
1065 */
1066real acosh(real x) @safe pure nothrow
1067{
1068    if (x > 1/real.epsilon)
1069        return LN2 + log(x);
1070    else
1071        return log(x + sqrt(x*x - 1));
1072}
1073
1074/// ditto
1075double acosh(double x) @safe pure nothrow { return acosh(cast(real)x); }
1076
1077/// ditto
1078float acosh(float x) @safe pure nothrow  { return acosh(cast(real)x); }
1079
1080
1081unittest
1082{
1083    assert(isNaN(acosh(0.9)));
1084    assert(isNaN(acosh(real.nan)));
1085    assert(acosh(1.0)==0.0);
1086    assert(acosh(real.infinity) == real.infinity);
1087    assert(isNaN(acosh(0.5)));
1088    assert(equalsDigit(acosh(cosh(3.0)), 3, useDigits));
1089}
1090
1091/***********************************
1092 * Calculates the inverse hyperbolic sine of x.
1093 *
1094 *  Mathematically,
1095 *  ---------------
1096 *  asinh(x) =  log( x + sqrt( x*x + 1 )) // if x >= +0
1097 *  asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
1098 *  -------------
1099 *
1100 *    $(TABLE_SV
1101 *    $(SVH x,                asinh(x)       )
1102 *    $(SV  $(NAN),           $(NAN)         )
1103 *    $(SV  $(PLUSMN)0,       $(PLUSMN)0      )
1104 *    $(SV  $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
1105 *    )
1106 */
1107real asinh(real x) @safe pure nothrow
1108{
1109    return (fabs(x) > 1 / real.epsilon)
1110       // beyond this point, x*x + 1 == x*x
1111       ?  copysign(LN2 + log(fabs(x)), x)
1112       // sqrt(x*x + 1) ==  1 + x * x / ( 1 + sqrt(x*x + 1) )
1113       : copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
1114}
1115
1116/// ditto
1117double asinh(double x) @safe pure nothrow { return asinh(cast(real)x); }
1118
1119/// ditto
1120float asinh(float x) @safe pure nothrow { return asinh(cast(real)x); }
1121
1122unittest
1123{
1124    assert(isIdentical(asinh(0.0), 0.0));
1125    assert(isIdentical(asinh(-0.0), -0.0));
1126    assert(asinh(real.infinity) == real.infinity);
1127    assert(asinh(-real.infinity) == -real.infinity);
1128    assert(isNaN(asinh(real.nan)));
1129    assert(equalsDigit(asinh(sinh(3.0)), 3, useDigits));
1130}
1131
1132/***********************************
1133 * Calculates the inverse hyperbolic tangent of x,
1134 * returning a value from ranging from -1 to 1.
1135 *
1136 * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
1137 *
1138 *
1139 * $(TABLE_DOMRG
1140 *  $(DOMAIN -$(INFIN)..$(INFIN))
1141 *  $(RANGE  -1..1) )
1142 * $(TABLE_SV
1143 *    $(SVH  x,     acosh(x) )
1144 *    $(SV  $(NAN), $(NAN) )
1145 *    $(SV  $(PLUSMN)0, $(PLUSMN)0)
1146 *    $(SV  -$(INFIN), -0)
1147 * )
1148 */
1149real atanh(real x) @safe pure nothrow
1150{
1151    // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
1152    return  0.5 * log1p( 2 * x / (1 - x) );
1153}
1154
1155/// ditto
1156double atanh(double x) @safe pure nothrow { return atanh(cast(real)x); }
1157
1158/// ditto
1159float atanh(float x) @safe pure nothrow { return atanh(cast(real)x); }
1160
1161
1162unittest
1163{
1164    assert(isIdentical(atanh(0.0), 0.0));
1165    assert(isIdentical(atanh(-0.0),-0.0));
1166    assert(isNaN(atanh(real.nan)));
1167    assert(isNaN(atanh(-real.infinity)));
1168    assert(atanh(0.0) == 0);
1169    assert(equalsDigit(atanh(tanh(0.5L)), 0.5, useDigits));
1170}
1171
1172/*****************************************
1173 * Returns x rounded to a long value using the current rounding mode.
1174 * If the integer value of x is
1175 * greater than long.max, the result is
1176 * indeterminate.
1177 */
1178long rndtol(real x) @safe pure nothrow;    /* intrinsic */
1179
1180
1181/*****************************************
1182 * Returns x rounded to a long value using the FE_TONEAREST rounding mode.
1183 * If the integer value of x is
1184 * greater than long.max, the result is
1185 * indeterminate.
1186 */
1187extern (C) real rndtonl(real x);
1188
1189/***************************************
1190 * Compute square root of x.
1191 *
1192 *      $(TABLE_SV
1193 *      $(TR $(TH x)         $(TH sqrt(x))   $(TH invalid?))
1194 *      $(TR $(TD -0.0)      $(TD -0.0)      $(TD no))
1195 *      $(TR $(TD $(LT)0.0)  $(TD $(NAN))    $(TD yes))
1196 *      $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
1197 *      )
1198 */
1199float sqrt(float x) @safe pure nothrow;    /* intrinsic */
1200
1201/// ditto
1202double sqrt(double x) @safe pure nothrow;  /* intrinsic */
1203
1204/// ditto
1205real sqrt(real x) @safe pure nothrow;      /* intrinsic */
1206
1207unittest
1208{
1209    //ctfe
1210    enum ZX80 = sqrt(7.0f);
1211    enum ZX81 = sqrt(7.0);
1212    enum ZX82 = sqrt(7.0L);
1213}
1214
1215creal sqrt(creal z) @safe pure nothrow
1216{
1217    creal c;
1218    real x,y,w,r;
1219
1220    if (z == 0)
1221    {
1222        c = 0 + 0i;
1223    }
1224    else
1225    {
1226        real z_re = z.re;
1227        real z_im = z.im;
1228
1229        x = fabs(z_re);
1230        y = fabs(z_im);
1231        if (x >= y)
1232        {
1233            r = y / x;
1234            w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
1235        }
1236        else
1237        {
1238            r = x / y;
1239            w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
1240        }
1241
1242        if (z_re >= 0)
1243        {
1244            c = w + (z_im / (w + w)) * 1.0i;
1245        }
1246        else
1247        {
1248            if (z_im < 0)
1249                w = -w;
1250            c = z_im / (w + w) + w * 1.0i;
1251        }
1252    }
1253    return c;
1254}
1255
1256/**
1257 * Calculates e$(SUP x).
1258 *
1259 *  $(TABLE_SV
1260 *    $(TR $(TH x)             $(TH e$(SUP x)) )
1261 *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
1262 *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
1263 *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
1264 *  )
1265 */
1266real exp(real x) @trusted pure nothrow
1267{
1268    version(D_InlineAsm_X86)
1269    {
1270        //  e^^x = 2^^(LOG2E*x)
1271        // (This is valid because the overflow & underflow limits for exp
1272        // and exp2 are so similar).
1273        return exp2(LOG2E*x);
1274    }
1275    else version(D_InlineAsm_X86_64)
1276    {
1277        //  e^^x = 2^^(LOG2E*x)
1278        // (This is valid because the overflow & underflow limits for exp
1279        // and exp2 are so similar).
1280        return exp2(LOG2E*x);
1281    }
1282    else
1283    {
1284        // Coefficients for exp(x)
1285        static immutable real[3] P = [
1286            9.9999999999999999991025E-1L,
1287            3.0299440770744196129956E-2L,
1288            1.2617719307481059087798E-4L,
1289        ];
1290        static immutable real[4] Q = [
1291            2.0000000000000000000897E0L,
1292            2.2726554820815502876593E-1L,
1293            2.5244834034968410419224E-3L,
1294            3.0019850513866445504159E-6L,
1295        ];
1296
1297        // C1 + C2 = LN2.
1298        enum real C1 = 6.9314575195312500000000E-1L;
1299        enum real C2 = 1.428606820309417232121458176568075500134E-6L;
1300
1301        // Overflow and Underflow limits.
1302        enum real OF =  1.1356523406294143949492E4L;
1303        enum real UF = -1.1432769596155737933527E4L;
1304
1305        // Special cases.
1306        if (isNaN(x))
1307            return x;
1308        if (x > OF)
1309            return real.infinity;
1310        if (x < UF)
1311            return 0.0;
1312
1313        // Express: e^^x = e^^g * 2^^n
1314        //   = e^^g * e^^(n * LOG2E)
1315        //   = e^^(g + n * LOG2E)
1316        int n = cast(int)floor(LOG2E * x + 0.5);
1317        x -= n * C1;
1318        x -= n * C2;
1319
1320        // Rational approximation for exponential of the fractional part:
1321        //  e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
1322        real xx = x * x;
1323        real px = x * poly(xx, P);
1324        x = px / (poly(xx, Q) - px);
1325        x = 1.0 + ldexp(x, 1);
1326
1327        // Scale by power of 2.
1328        x = ldexp(x, n);
1329
1330        return x;
1331    }
1332}
1333
1334/// ditto
1335double exp(double x) @safe pure nothrow  { return exp(cast(real)x); }
1336
1337/// ditto
1338float exp(float x)  @safe pure nothrow   { return exp(cast(real)x); }
1339
1340unittest
1341{
1342    assert(equalsDigit(exp(3.0L), E * E * E, useDigits));
1343}
1344
1345/**
1346 * Calculates the value of the natural logarithm base (e)
1347 * raised to the power of x, minus 1.
1348 *
1349 * For very small x, expm1(x) is more accurate
1350 * than exp(x)-1.
1351 *
1352 *  $(TABLE_SV
1353 *    $(TR $(TH x)             $(TH e$(SUP x)-1)  )
1354 *    $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) )
1355 *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN))    )
1356 *    $(TR $(TD -$(INFIN))     $(TD -1.0)         )
1357 *    $(TR $(TD $(NAN))        $(TD $(NAN))       )
1358 *  )
1359 */
1360real expm1(real x) @trusted pure nothrow
1361{
1362    version(D_InlineAsm_X86)
1363    {
1364        enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1365        asm
1366        {
1367            /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1368             * Author: Don Clugston.
1369             *
1370             *    expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
1371             *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
1372             *     and 2ym1 = (2^^(y-rndint(y))-1).
1373             *    If 2rndy  < 0.5*real.epsilon, result is -1.
1374             *    Implementation is otherwise the same as for exp2()
1375             */
1376            naked;
1377            fld real ptr [ESP+4] ; // x
1378            mov AX, [ESP+4+8]; // AX = exponent and sign
1379            sub ESP, 12+8; // Create scratch space on the stack
1380            // [ESP,ESP+2] = scratchint
1381            // [ESP+4..+6, +8..+10, +10] = scratchreal
1382            // set scratchreal mantissa = 1.0
1383            mov dword ptr [ESP+8], 0;
1384            mov dword ptr [ESP+8+4], 0x80000000;
1385            and AX, 0x7FFF; // drop sign bit
1386            cmp AX, 0x401D; // avoid InvalidException in fist
1387            jae L_extreme;
1388            fldl2e;
1389            fmulp ST(1), ST; // y = x*log2(e)
1390            fist dword ptr [ESP]; // scratchint = rndint(y)
1391            fisub dword ptr [ESP]; // y - rndint(y)
1392            // and now set scratchreal exponent
1393            mov EAX, [ESP];
1394            add EAX, 0x3fff;
1395            jle short L_largenegative;
1396            cmp EAX,0x8000;
1397            jge short L_largepositive;
1398            mov [ESP+8+8],AX;
1399            f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
1400            fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
1401            fmul ST(1), ST;  // ST=2rndy, ST(1)=2rndy*2ym1
1402            fld1;
1403            fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
1404            faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
1405            add ESP,12+8;
1406            ret PARAMSIZE;
1407
1408L_extreme:  // Extreme exponent. X is very large positive, very
1409            // large negative, infinity, or NaN.
1410            fxam;
1411            fstsw AX;
1412            test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1413            jz L_was_nan;  // if x is NaN, returns x
1414            test AX, 0x0200;
1415            jnz L_largenegative;
1416L_largepositive:
1417            // Set scratchreal = real.max.
1418            // squaring it will create infinity, and set overflow flag.
1419            mov word  ptr [ESP+8+8], 0x7FFE;
1420            fstp ST(0);
1421            fld real ptr [ESP+8];  // load scratchreal
1422            fmul ST(0), ST;        // square it, to create havoc!
1423L_was_nan:
1424            add ESP,12+8;
1425            ret PARAMSIZE;
1426L_largenegative:
1427            fstp ST(0);
1428            fld1;
1429            fchs; // return -1. Underflow flag is not set.
1430            add ESP,12+8;
1431            ret PARAMSIZE;
1432        }
1433    }
1434    else version(D_InlineAsm_X86_64)
1435    {
1436        asm
1437        {
1438            naked;
1439        }
1440        version (Win64)
1441        {
1442            asm
1443            {
1444                fld   real ptr [RCX];  // x
1445                mov   AX,[RCX+8];      // AX = exponent and sign
1446            }
1447        }
1448        else
1449        {
1450            asm
1451            {
1452                fld   real ptr [RSP+8];  // x
1453                mov   AX,[RSP+8+8];      // AX = exponent and sign
1454            }
1455        }
1456        asm
1457        {
1458            /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1459             * Author: Don Clugston.
1460             *
1461             *    expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
1462             *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
1463             *     and 2ym1 = (2^(y-rndint(y))-1).
1464             *    If 2rndy  < 0.5*real.epsilon, result is -1.
1465             *    Implementation is otherwise the same as for exp2()
1466             */
1467            sub RSP, 24;       // Create scratch space on the stack
1468            // [RSP,RSP+2] = scratchint
1469            // [RSP+4..+6, +8..+10, +10] = scratchreal
1470            // set scratchreal mantissa = 1.0
1471            mov dword ptr [RSP+8], 0;
1472            mov dword ptr [RSP+8+4], 0x80000000;
1473            and AX, 0x7FFF; // drop sign bit
1474            cmp AX, 0x401D; // avoid InvalidException in fist
1475            jae L_extreme;
1476            fldl2e;
1477            fmul ; // y = x*log2(e)
1478            fist dword ptr [RSP]; // scratchint = rndint(y)
1479            fisub dword ptr [RSP]; // y - rndint(y)
1480            // and now set scratchreal exponent
1481            mov EAX, [RSP];
1482            add EAX, 0x3fff;
1483            jle short L_largenegative;
1484            cmp EAX,0x8000;
1485            jge short L_largepositive;
1486            mov [RSP+8+8],AX;
1487            f2xm1; // 2^(y-rndint(y)) -1
1488            fld real ptr [RSP+8] ; // 2^rndint(y)
1489            fmul ST(1), ST;
1490            fld1;
1491            fsubp ST(1), ST;
1492            fadd;
1493            add RSP,24;
1494            ret;
1495
1496L_extreme:  // Extreme exponent. X is very large positive, very
1497            // large negative, infinity, or NaN.
1498            fxam;
1499            fstsw AX;
1500            test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1501            jz L_was_nan;  // if x is NaN, returns x
1502            test AX, 0x0200;
1503            jnz L_largenegative;
1504L_largepositive:
1505            // Set scratchreal = real.max.
1506            // squaring it will create infinity, and set overflow flag.
1507            mov word  ptr [RSP+8+8], 0x7FFE;
1508            fstp ST(0);
1509            fld real ptr [RSP+8];  // load scratchreal
1510            fmul ST(0), ST;        // square it, to create havoc!
1511L_was_nan:
1512            add RSP,24;
1513            ret;
1514
1515L_largenegative:
1516            fstp ST(0);
1517            fld1;
1518            fchs; // return -1. Underflow flag is not set.
1519            add RSP,24;
1520            ret;
1521        }
1522    }
1523    else
1524    {
1525        // Coefficients for exp(x) - 1
1526        static immutable real[5] P = [
1527           -1.586135578666346600772998894928250240826E4L,
1528            2.642771505685952966904660652518429479531E3L,
1529           -3.423199068835684263987132888286791620673E2L,
1530            1.800826371455042224581246202420972737840E1L,
1531           -5.238523121205561042771939008061958820811E-1L,
1532        ];
1533        static immutable real[6] Q = [
1534           -9.516813471998079611319047060563358064497E4L,
1535            3.964866271411091674556850458227710004570E4L,
1536           -7.207678383830091850230366618190187434796E3L,
1537            7.206038318724600171970199625081491823079E2L,
1538           -4.002027679107076077238836622982900945173E1L,
1539            1.000000000000000000000000000000000000000E0L,
1540        ];
1541
1542        // C1 + C2 = LN2.
1543        enum real C1 = 6.9314575195312500000000E-1L;
1544        enum real C2 = 1.4286068203094172321215E-6L;
1545
1546        // Overflow and Underflow limits.
1547        enum real OF =  1.1356523406294143949492E4L;
1548        enum real UF = -4.5054566736396445112120088E1L;
1549
1550        // Special cases.
1551        if (x > OF)
1552            return real.infinity;
1553        if (x == 0.0)
1554            return x;
1555        if (x < UF)
1556            return -1.0;
1557
1558        // Express x = LN2 (n + remainder), remainder not exceeding 1/2.
1559        int n = cast(int)floor(0.5 + x / LN2);
1560        x -= n * C1;
1561        x -= n * C2;
1562
1563        // Rational approximation:
1564        //  exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x)
1565        real px = x * poly(x, P);
1566        real qx = poly(x, Q);
1567        real xx = x * x;
1568        qx = x + (0.5 * xx + xx * px / qx);
1569
1570        // We have qx = exp(remainder LN2) - 1, so:
1571        //  exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1.
1572        px = ldexp(1.0, n);
1573        x = px * qx + (px - 1.0);
1574
1575        return x;
1576    }
1577}
1578
1579
1580
1581/**
1582 * Calculates 2$(SUP x).
1583 *
1584 *  $(TABLE_SV
1585 *    $(TR $(TH x)             $(TH exp2(x))   )
1586 *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
1587 *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
1588 *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
1589 *  )
1590 */
1591real exp2(real x) @trusted pure nothrow
1592{
1593    version(D_InlineAsm_X86)
1594    {
1595        enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1596
1597        asm
1598        {
1599            /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1600             * Author: Don Clugston.
1601             *
1602             * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
1603             * The trick for high performance is to avoid the fscale(28cycles on core2),
1604             * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1605             *
1606             * We can do frndint by using fist. BUT we can't use it for huge numbers,
1607             * because it will set the Invalid Operation flag if overflow or NaN occurs.
1608             * Fortunately, whenever this happens the result would be zero or infinity.
1609             *
1610             * We can perform fscale by directly poking into the exponent. BUT this doesn't
1611             * work for the (very rare) cases where the result is subnormal. So we fall back
1612             * to the slow method in that case.
1613             */
1614            naked;
1615            fld real ptr [ESP+4] ; // x
1616            mov AX, [ESP+4+8]; // AX = exponent and sign
1617            sub ESP, 12+8; // Create scratch space on the stack
1618            // [ESP,ESP+2] = scratchint
1619            // [ESP+4..+6, +8..+10, +10] = scratchreal
1620            // set scratchreal mantissa = 1.0
1621            mov dword ptr [ESP+8], 0;
1622            mov dword ptr [ESP+8+4], 0x80000000;
1623            and AX, 0x7FFF; // drop sign bit
1624            cmp AX, 0x401D; // avoid InvalidException in fist
1625            jae L_extreme;
1626            fist dword ptr [ESP]; // scratchint = rndint(x)
1627            fisub dword ptr [ESP]; // x - rndint(x)
1628            // and now set scratchreal exponent
1629            mov EAX, [ESP];
1630            add EAX, 0x3fff;
1631            jle short L_subnormal;
1632            cmp EAX,0x8000;
1633            jge short L_overflow;
1634            mov [ESP+8+8],AX;
1635L_normal:
1636            f2xm1;
1637            fld1;
1638            faddp ST(1), ST; // 2^^(x-rndint(x))
1639            fld real ptr [ESP+8] ; // 2^^rndint(x)
1640            add ESP,12+8;
1641            fmulp ST(1), ST;
1642            ret PARAMSIZE;
1643
1644L_subnormal:
1645            // Result will be subnormal.
1646            // In this rare case, the simple poking method doesn't work.
1647            // The speed doesn't matter, so use the slow fscale method.
1648            fild dword ptr [ESP];  // scratchint
1649            fld1;
1650            fscale;
1651            fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
1652            fstp ST(0);         // drop scratchint
1653            jmp L_normal;
1654
1655L_extreme:  // Extreme exponent. X is very large positive, very
1656            // large negative, infinity, or NaN.
1657            fxam;
1658            fstsw AX;
1659            test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1660            jz L_was_nan;  // if x is NaN, returns x
1661            // set scratchreal = real.min_normal
1662            // squaring it will return 0, setting underflow flag
1663            mov word  ptr [ESP+8+8], 1;
1664            test AX, 0x0200;
1665            jnz L_waslargenegative;
1666L_overflow:
1667            // Set scratchreal = real.max.
1668            // squaring it will create infinity, and set overflow flag.
1669            mov word  ptr [ESP+8+8], 0x7FFE;
1670L_waslargenegative:
1671            fstp ST(0);
1672            fld real ptr [ESP+8];  // load scratchreal
1673            fmul ST(0), ST;        // square it, to create havoc!
1674L_was_nan:
1675            add ESP,12+8;
1676            ret PARAMSIZE;
1677        }
1678    }
1679    else version(D_InlineAsm_X86_64)
1680    {
1681        asm
1682        {
1683            naked;
1684        }
1685        version (Win64)
1686        {
1687            asm
1688            {
1689                fld   real ptr [RCX];  // x
1690                mov   AX,[RCX+8];      // AX = exponent and sign
1691            }
1692        }
1693        else
1694        {
1695            asm
1696            {
1697                fld   real ptr [RSP+8];  // x
1698                mov   AX,[RSP+8+8];      // AX = exponent and sign
1699            }
1700        }
1701        asm
1702        {
1703            /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1704             * Author: Don Clugston.
1705             *
1706             * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
1707             * The trick for high performance is to avoid the fscale(28cycles on core2),
1708             * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1709             *
1710             * We can do frndint by using fist. BUT we can't use it for huge numbers,
1711             * because it will set the Invalid Operation flag is overflow or NaN occurs.
1712             * Fortunately, whenever this happens the result would be zero or infinity.
1713             *
1714             * We can perform fscale by directly poking into the exponent. BUT this doesn't
1715             * work f

Large files files are truncated, but you can click here to view the full file