/std/math.d
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 = Γ 33 * THETA = θ 34 * INTEGRAL = ∫ 35 * INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>) 36 * POWER = $1<sup>$2</sup> 37 * SUB = $1<sub>$2</sub> 38 * BIGSUM = $(BIG Σ <sup>$2</sup><sub>$(SMALL $1)</sub>) 39 * CHOOSE = $(BIG () <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG )) 40 * PLUSMN = ± 41 * INFIN = ∞ 42 * PLUSMNINF = ±∞ 43 * PI = π 44 * LT = < 45 * GT = > 46 * SQRT = √ 47 * HALF = ½ 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