1 // Written in the D programming language. 2 3 /** 4 This is a submodule of $(MREF std, math). 5 6 It contains several exponential and logarithm functions. 7 8 Copyright: Copyright The D Language Foundation 2000 - 2011. 9 D implementations of exp, expm1, exp2, log, log10, log1p, and log2 10 functions are based on the CEPHES math library, which is Copyright 11 (C) 2001 Stephen L. Moshier $(LT)steve (at) moshier.net$(GT) and are 12 incorporated herein by permission of the author. The author reserves 13 the right to distribute this material elsewhere under different 14 copying permissions. These modifications are distributed here under 15 the following terms: 16 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). 17 Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston, 18 Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger 19 Source: $(PHOBOSSRC std/math/exponential.d) 20 21 Macros: 22 TABLE_SV = <table border="1" cellpadding="4" cellspacing="0"> 23 <caption>Special Values</caption> 24 $0</table> 25 NAN = $(RED NAN) 26 PLUSMN = ± 27 INFIN = ∞ 28 PLUSMNINF = ±∞ 29 LT = < 30 GT = > 31 */ 32 33 module std.math.exponential; 34 35 import std.traits : isFloatingPoint, isIntegral, isSigned, isUnsigned, Largest, Unqual; 36 37 static import core.math; 38 static import core.stdc.math; 39 40 version (DigitalMars) 41 { 42 version = INLINE_YL2X; // x87 has opcodes for these 43 } 44 45 version (D_InlineAsm_X86) version = InlineAsm_X86_Any; 46 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any; 47 48 version (InlineAsm_X86_Any) version = InlineAsm_X87; 49 version (InlineAsm_X87) 50 { 51 static assert(real.mant_dig == 64); 52 version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC; 53 } 54 55 version (D_HardFloat) 56 { 57 // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport 58 version (IeeeFlagsSupport) version = FloatingPointControlSupport; 59 } 60 61 /** 62 * Compute the value of x $(SUPERSCRIPT n), where n is an integer 63 */ 64 Unqual!F pow(F, G)(F x, G n) @nogc @trusted pure nothrow 65 if (isFloatingPoint!(F) && isIntegral!(G)) 66 { 67 import std.traits : Unsigned; 68 69 real p = 1.0, v = void; 70 Unsigned!(Unqual!G) m = n; 71 72 if (n < 0) 73 { 74 if (n == -1) return 1 / x; 75 76 m = cast(typeof(m))(0 - n); 77 v = p / x; 78 } 79 else 80 { 81 switch (n) 82 { 83 case 0: 84 return 1.0; 85 case 1: 86 return x; 87 case 2: 88 return x * x; 89 default: 90 } 91 92 v = x; 93 } 94 95 while (1) 96 { 97 if (m & 1) 98 p *= v; 99 m >>= 1; 100 if (!m) 101 break; 102 v *= v; 103 } 104 return p; 105 } 106 107 /// 108 @safe pure nothrow @nogc unittest 109 { 110 import std.math.operations : feqrel; 111 112 assert(pow(2.0, 5) == 32.0); 113 assert(pow(1.5, 9).feqrel(38.4433) > 16); 114 assert(pow(real.nan, 2) is real.nan); 115 assert(pow(real.infinity, 2) == real.infinity); 116 } 117 118 @safe pure nothrow @nogc unittest 119 { 120 import std.math.operations : isClose, feqrel; 121 122 // Make sure it instantiates and works properly on immutable values and 123 // with various integer and float types. 124 immutable real x = 46; 125 immutable float xf = x; 126 immutable double xd = x; 127 immutable uint one = 1; 128 immutable ushort two = 2; 129 immutable ubyte three = 3; 130 immutable ulong eight = 8; 131 132 immutable int neg1 = -1; 133 immutable short neg2 = -2; 134 immutable byte neg3 = -3; 135 immutable long neg8 = -8; 136 137 138 assert(pow(x,0) == 1.0); 139 assert(pow(xd,one) == x); 140 assert(pow(xf,two) == x * x); 141 assert(pow(x,three) == x * x * x); 142 assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x)); 143 144 assert(pow(x, neg1) == 1 / x); 145 146 assert(isClose(pow(xd, neg2), cast(double) (1 / (x * x)), 1e-15)); 147 assert(isClose(pow(xf, neg8), cast(float) (1 / ((x * x) * (x * x) * (x * x) * (x * x))), 1e-15)); 148 149 assert(feqrel(pow(x, neg3), 1 / (x * x * x)) >= real.mant_dig - 1); 150 } 151 152 @safe @nogc nothrow unittest 153 { 154 import std.math.operations : isClose; 155 156 assert(isClose(pow(2.0L, 10L), 1024, 1e-18)); 157 } 158 159 // https://issues.dlang.org/show_bug.cgi?id=21601 160 @safe @nogc nothrow pure unittest 161 { 162 // When reals are large enough the results of pow(b, e) can be 163 // calculated correctly, if b is of type float or double and e is 164 // not too large. 165 static if (real.mant_dig >= 64) 166 { 167 // expected result: 3.790e-42 168 assert(pow(-513645318757045764096.0f, -2) > 0.0); 169 170 // expected result: 3.763915357831797e-309 171 assert(pow(-1.6299717435255677e+154, -2) > 0.0); 172 } 173 } 174 175 @safe @nogc nothrow unittest 176 { 177 import std.math.operations : isClose; 178 import std.math.traits : isInfinity; 179 180 static float f1 = 19100.0f; 181 static float f2 = 0.000012f; 182 183 assert(isClose(pow(f1,9), 3.3829868e+38f)); 184 assert(isInfinity(pow(f1,10))); 185 assert(pow(f2,9) > 0.0f); 186 assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal)); 187 188 static double d1 = 21800.0; 189 static double d2 = 0.000012; 190 191 assert(isClose(pow(d1,71), 1.0725339442974e+308)); 192 assert(isInfinity(pow(d1,72))); 193 assert(pow(d2,65) > 0.0f); 194 assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal)); 195 196 static if (real.mant_dig == 64) // x87 197 { 198 static real r1 = 21950.0L; 199 static real r2 = 0.000011L; 200 201 assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L)); 202 assert(isInfinity(pow(r1,1137))); 203 assert(pow(r2,998) > 0.0L); 204 assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal)); 205 } 206 } 207 208 @safe @nogc nothrow pure unittest 209 { 210 import std.math.operations : isClose; 211 212 enum f1 = 19100.0f; 213 enum f2 = 0.000012f; 214 215 static assert(isClose(pow(f1,9), 3.3829868e+38f)); 216 static assert(pow(f1,10) > float.max); 217 static assert(pow(f2,9) > 0.0f); 218 static assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal)); 219 220 enum d1 = 21800.0; 221 enum d2 = 0.000012; 222 223 static assert(isClose(pow(d1,71), 1.0725339442974e+308)); 224 static assert(pow(d1,72) > double.max); 225 static assert(pow(d2,65) > 0.0f); 226 static assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal)); 227 228 static if (real.mant_dig == 64) // x87 229 { 230 enum r1 = 21950.0L; 231 enum r2 = 0.000011L; 232 233 static assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L)); 234 static assert(pow(r1,1137) > real.max); 235 static assert(pow(r2,998) > 0.0L); 236 static assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal)); 237 } 238 } 239 240 /** 241 * Compute the power of two integral numbers. 242 * 243 * Params: 244 * x = base 245 * n = exponent 246 * 247 * Returns: 248 * x raised to the power of n. If n is negative the result is 1 / pow(x, -n), 249 * which is calculated as integer division with remainder. This may result in 250 * a division by zero error. 251 * 252 * If both x and n are 0, the result is 1. 253 * 254 * Throws: 255 * If x is 0 and n is negative, the result is the same as the result of a 256 * division by zero. 257 */ 258 typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @nogc @trusted pure nothrow 259 if (isIntegral!(F) && isIntegral!(G)) 260 { 261 import std.traits : isSigned; 262 263 typeof(return) p, v = void; 264 Unqual!G m = n; 265 266 static if (isSigned!(F)) 267 { 268 if (x == -1) return cast(typeof(return)) (m & 1 ? -1 : 1); 269 } 270 static if (isSigned!(G)) 271 { 272 if (x == 0 && m <= -1) return x / 0; 273 } 274 if (x == 1) return 1; 275 static if (isSigned!(G)) 276 { 277 if (m < 0) return 0; 278 } 279 280 switch (m) 281 { 282 case 0: 283 p = 1; 284 break; 285 286 case 1: 287 p = x; 288 break; 289 290 case 2: 291 p = x * x; 292 break; 293 294 default: 295 v = x; 296 p = 1; 297 while (1) 298 { 299 if (m & 1) 300 p *= v; 301 m >>= 1; 302 if (!m) 303 break; 304 v *= v; 305 } 306 break; 307 } 308 return p; 309 } 310 311 /// 312 @safe pure nothrow @nogc unittest 313 { 314 assert(pow(2, 3) == 8); 315 assert(pow(3, 2) == 9); 316 317 assert(pow(2, 10) == 1_024); 318 assert(pow(2, 20) == 1_048_576); 319 assert(pow(2, 30) == 1_073_741_824); 320 321 assert(pow(0, 0) == 1); 322 323 assert(pow(1, -5) == 1); 324 assert(pow(1, -6) == 1); 325 assert(pow(-1, -5) == -1); 326 assert(pow(-1, -6) == 1); 327 328 assert(pow(-2, 5) == -32); 329 assert(pow(-2, -5) == 0); 330 assert(pow(cast(double) -2, -5) == -0.03125); 331 } 332 333 @safe pure nothrow @nogc unittest 334 { 335 immutable int one = 1; 336 immutable byte two = 2; 337 immutable ubyte three = 3; 338 immutable short four = 4; 339 immutable long ten = 10; 340 341 assert(pow(two, three) == 8); 342 assert(pow(two, ten) == 1024); 343 assert(pow(one, ten) == 1); 344 assert(pow(ten, four) == 10_000); 345 assert(pow(four, 10) == 1_048_576); 346 assert(pow(three, four) == 81); 347 } 348 349 // https://issues.dlang.org/show_bug.cgi?id=7006 350 @safe pure nothrow @nogc unittest 351 { 352 assert(pow(5, -1) == 0); 353 assert(pow(-5, -1) == 0); 354 assert(pow(5, -2) == 0); 355 assert(pow(-5, -2) == 0); 356 assert(pow(-1, int.min) == 1); 357 assert(pow(-2, int.min) == 0); 358 359 assert(pow(4294967290UL,2) == 18446744022169944100UL); 360 assert(pow(0,uint.max) == 0); 361 } 362 363 /**Computes integer to floating point powers.*/ 364 real pow(I, F)(I x, F y) @nogc @trusted pure nothrow 365 if (isIntegral!I && isFloatingPoint!F) 366 { 367 return pow(cast(real) x, cast(Unqual!F) y); 368 } 369 370 /// 371 @safe pure nothrow @nogc unittest 372 { 373 assert(pow(2, 5.0) == 32.0); 374 assert(pow(7, 3.0) == 343.0); 375 assert(pow(2, real.nan) is real.nan); 376 assert(pow(2, real.infinity) == real.infinity); 377 } 378 379 /** 380 * Calculates x$(SUPERSCRIPT y). 381 * 382 * $(TABLE_SV 383 * $(TR $(TH x) $(TH y) $(TH pow(x, y)) 384 * $(TH div 0) $(TH invalid?)) 385 * $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD 1.0) 386 * $(TD no) $(TD no) ) 387 * $(TR $(TD |x| $(GT) 1) $(TD +$(INFIN)) $(TD +$(INFIN)) 388 * $(TD no) $(TD no) ) 389 * $(TR $(TD |x| $(LT) 1) $(TD +$(INFIN)) $(TD +0.0) 390 * $(TD no) $(TD no) ) 391 * $(TR $(TD |x| $(GT) 1) $(TD -$(INFIN)) $(TD +0.0) 392 * $(TD no) $(TD no) ) 393 * $(TR $(TD |x| $(LT) 1) $(TD -$(INFIN)) $(TD +$(INFIN)) 394 * $(TD no) $(TD no) ) 395 * $(TR $(TD +$(INFIN)) $(TD $(GT) 0.0) $(TD +$(INFIN)) 396 * $(TD no) $(TD no) ) 397 * $(TR $(TD +$(INFIN)) $(TD $(LT) 0.0) $(TD +0.0) 398 * $(TD no) $(TD no) ) 399 * $(TR $(TD -$(INFIN)) $(TD odd integer $(GT) 0.0) $(TD -$(INFIN)) 400 * $(TD no) $(TD no) ) 401 * $(TR $(TD -$(INFIN)) $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN)) 402 * $(TD no) $(TD no)) 403 * $(TR $(TD -$(INFIN)) $(TD odd integer $(LT) 0.0) $(TD -0.0) 404 * $(TD no) $(TD no) ) 405 * $(TR $(TD -$(INFIN)) $(TD $(LT) 0.0, not odd integer) $(TD +0.0) 406 * $(TD no) $(TD no) ) 407 * $(TR $(TD $(PLUSMN)1.0) $(TD $(PLUSMN)$(INFIN)) $(TD -$(NAN)) 408 * $(TD no) $(TD yes) ) 409 * $(TR $(TD $(LT) 0.0) $(TD finite, nonintegral) $(TD $(NAN)) 410 * $(TD no) $(TD yes)) 411 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(LT) 0.0) $(TD $(PLUSMNINF)) 412 * $(TD yes) $(TD no) ) 413 * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN)) 414 * $(TD yes) $(TD no)) 415 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(GT) 0.0) $(TD $(PLUSMN)0.0) 416 * $(TD no) $(TD no) ) 417 * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT) 0.0, not odd integer) $(TD +0.0) 418 * $(TD no) $(TD no) ) 419 * ) 420 */ 421 Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow 422 if (isFloatingPoint!(F) && isFloatingPoint!(G)) 423 { 424 import core.math : fabs, sqrt; 425 import std.math.traits : isInfinity, isNaN, signbit; 426 427 alias Float = typeof(return); 428 429 static real impl(real x, real y) @nogc pure nothrow 430 { 431 // Special cases. 432 if (isNaN(y)) 433 return y; 434 if (isNaN(x) && y != 0.0) 435 return x; 436 437 // Even if x is NaN. 438 if (y == 0.0) 439 return 1.0; 440 if (y == 1.0) 441 return x; 442 443 if (isInfinity(y)) 444 { 445 if (isInfinity(x)) 446 { 447 if (!signbit(y) && !signbit(x)) 448 return F.infinity; 449 else 450 return F.nan; 451 } 452 else if (fabs(x) > 1) 453 { 454 if (signbit(y)) 455 return +0.0; 456 else 457 return F.infinity; 458 } 459 else if (fabs(x) == 1) 460 { 461 return F.nan; 462 } 463 else // < 1 464 { 465 if (signbit(y)) 466 return F.infinity; 467 else 468 return +0.0; 469 } 470 } 471 if (isInfinity(x)) 472 { 473 if (signbit(x)) 474 { 475 long i = cast(long) y; 476 if (y > 0.0) 477 { 478 if (i == y && i & 1) 479 return -F.infinity; 480 else if (i == y) 481 return F.infinity; 482 else 483 return -F.nan; 484 } 485 else if (y < 0.0) 486 { 487 if (i == y && i & 1) 488 return -0.0; 489 else if (i == y) 490 return +0.0; 491 else 492 return F.nan; 493 } 494 } 495 else 496 { 497 if (y > 0.0) 498 return F.infinity; 499 else if (y < 0.0) 500 return +0.0; 501 } 502 } 503 504 if (x == 0.0) 505 { 506 if (signbit(x)) 507 { 508 long i = cast(long) y; 509 if (y > 0.0) 510 { 511 if (i == y && i & 1) 512 return -0.0; 513 else 514 return +0.0; 515 } 516 else if (y < 0.0) 517 { 518 if (i == y && i & 1) 519 return -F.infinity; 520 else 521 return F.infinity; 522 } 523 } 524 else 525 { 526 if (y > 0.0) 527 return +0.0; 528 else if (y < 0.0) 529 return F.infinity; 530 } 531 } 532 if (x == 1.0) 533 return 1.0; 534 535 if (y >= F.max) 536 { 537 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0)) 538 return 0.0; 539 if (x > 1.0 || x < -1.0) 540 return F.infinity; 541 } 542 if (y <= -F.max) 543 { 544 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0)) 545 return F.infinity; 546 if (x > 1.0 || x < -1.0) 547 return 0.0; 548 } 549 550 if (x >= F.max) 551 { 552 if (y > 0.0) 553 return F.infinity; 554 else 555 return 0.0; 556 } 557 if (x <= -F.max) 558 { 559 long i = cast(long) y; 560 if (y > 0.0) 561 { 562 if (i == y && i & 1) 563 return -F.infinity; 564 else 565 return F.infinity; 566 } 567 else if (y < 0.0) 568 { 569 if (i == y && i & 1) 570 return -0.0; 571 else 572 return +0.0; 573 } 574 } 575 576 // Integer power of x. 577 long iy = cast(long) y; 578 if (iy == y && fabs(y) < 32_768.0) 579 return pow(x, iy); 580 581 real sign = 1.0; 582 if (x < 0) 583 { 584 // Result is real only if y is an integer 585 // Check for a non-zero fractional part 586 enum maxOdd = pow(2.0L, real.mant_dig) - 1.0L; 587 static if (maxOdd > ulong.max) 588 { 589 // Generic method, for any FP type 590 import std.math.rounding : floor; 591 if (floor(y) != y) 592 return sqrt(x); // Complex result -- create a NaN 593 594 const hy = 0.5 * y; 595 if (floor(hy) != hy) 596 sign = -1.0; 597 } 598 else 599 { 600 // Much faster, if ulong has enough precision 601 const absY = fabs(y); 602 if (absY <= maxOdd) 603 { 604 const uy = cast(ulong) absY; 605 if (uy != absY) 606 return sqrt(x); // Complex result -- create a NaN 607 608 if (uy & 1) 609 sign = -1.0; 610 } 611 } 612 x = -x; 613 } 614 version (INLINE_YL2X) 615 { 616 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) ) 617 // TODO: This is not accurate in practice. A fast and accurate 618 // (though complicated) method is described in: 619 // "An efficient rounding boundary test for pow(x, y) 620 // in double precision", C.Q. Lauter and V. Lefvre, INRIA (2007). 621 return sign * exp2( core.math.yl2x(x, y) ); 622 } 623 else 624 { 625 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) ) 626 // TODO: This is not accurate in practice. A fast and accurate 627 // (though complicated) method is described in: 628 // "An efficient rounding boundary test for pow(x, y) 629 // in double precision", C.Q. Lauter and V. Lefvre, INRIA (2007). 630 Float w = exp2(y * log2(x)); 631 return sign * w; 632 } 633 } 634 return impl(x, y); 635 } 636 637 /// 638 @safe pure nothrow @nogc unittest 639 { 640 import std.math.operations : isClose; 641 642 assert(isClose(pow(2.0, 3.0), 8.0)); 643 assert(isClose(pow(1.5, 10.0), 57.6650390625)); 644 645 // square root of 9 646 assert(isClose(pow(9.0, 0.5), 3.0)); 647 // 10th root of 1024 648 assert(isClose(pow(1024.0, 0.1), 2.0)); 649 650 assert(isClose(pow(-4.0, 3.0), -64.0)); 651 652 // reciprocal of 4 ^^ 2 653 assert(isClose(pow(4.0, -2.0), 0.0625)); 654 // reciprocal of (-2) ^^ 3 655 assert(isClose(pow(-2.0, -3.0), -0.125)); 656 657 assert(isClose(pow(-2.5, 3.0), -15.625)); 658 // reciprocal of 2.5 ^^ 3 659 assert(isClose(pow(2.5, -3.0), 0.064)); 660 // reciprocal of (-2.5) ^^ 3 661 assert(isClose(pow(-2.5, -3.0), -0.064)); 662 663 // reciprocal of square root of 4 664 assert(isClose(pow(4.0, -0.5), 0.5)); 665 666 // per definition 667 assert(isClose(pow(0.0, 0.0), 1.0)); 668 } 669 670 /// 671 @safe pure nothrow @nogc unittest 672 { 673 import std.math.operations : isClose; 674 675 // the result is a complex number 676 // which cannot be represented as floating point number 677 import std.math.traits : isNaN; 678 assert(isNaN(pow(-2.5, -1.5))); 679 680 // use the ^^-operator of std.complex instead 681 import std.complex : complex; 682 auto c1 = complex(-2.5, 0.0); 683 auto c2 = complex(-1.5, 0.0); 684 auto result = c1 ^^ c2; 685 // exact result apparently depends on `real` precision => increased tolerance 686 assert(isClose(result.re, -4.64705438e-17, 2e-4)); 687 assert(isClose(result.im, 2.52982e-1, 2e-4)); 688 } 689 690 @safe pure nothrow @nogc unittest 691 { 692 import std.math.traits : isNaN; 693 694 assert(pow(1.5, real.infinity) == real.infinity); 695 assert(pow(0.5, real.infinity) == 0.0); 696 assert(pow(1.5, -real.infinity) == 0.0); 697 assert(pow(0.5, -real.infinity) == real.infinity); 698 assert(pow(real.infinity, 1.0) == real.infinity); 699 assert(pow(real.infinity, -1.0) == 0.0); 700 assert(pow(real.infinity, real.infinity) == real.infinity); 701 assert(pow(-real.infinity, 1.0) == -real.infinity); 702 assert(pow(-real.infinity, 2.0) == real.infinity); 703 assert(pow(-real.infinity, -1.0) == -0.0); 704 assert(pow(-real.infinity, -2.0) == 0.0); 705 assert(isNaN(pow(1.0, real.infinity))); 706 assert(pow(0.0, -1.0) == real.infinity); 707 assert(pow(real.nan, 0.0) == 1.0); 708 assert(isNaN(pow(real.nan, 3.0))); 709 assert(isNaN(pow(3.0, real.nan))); 710 } 711 712 @safe @nogc nothrow unittest 713 { 714 import std.math.operations : isClose; 715 716 assert(isClose(pow(2.0L, 10.0L), 1024, 1e-18)); 717 } 718 719 @safe pure nothrow @nogc unittest 720 { 721 import std.math.operations : isClose; 722 import std.math.traits : isIdentical, isNaN; 723 import std.math.constants : PI; 724 725 // Test all the special values. These unittests can be run on Windows 726 // by temporarily changing the version (linux) to version (all). 727 immutable float zero = 0; 728 immutable real one = 1; 729 immutable double two = 2; 730 immutable float three = 3; 731 immutable float fnan = float.nan; 732 immutable double dnan = double.nan; 733 immutable real rnan = real.nan; 734 immutable dinf = double.infinity; 735 immutable rninf = -real.infinity; 736 737 assert(pow(fnan, zero) == 1); 738 assert(pow(dnan, zero) == 1); 739 assert(pow(rnan, zero) == 1); 740 741 assert(pow(two, dinf) == double.infinity); 742 assert(isIdentical(pow(0.2f, dinf), +0.0)); 743 assert(pow(0.99999999L, rninf) == real.infinity); 744 assert(isIdentical(pow(1.000000001, rninf), +0.0)); 745 assert(pow(dinf, 0.001) == dinf); 746 assert(isIdentical(pow(dinf, -0.001), +0.0)); 747 assert(pow(rninf, 3.0L) == rninf); 748 assert(pow(rninf, 2.0L) == real.infinity); 749 assert(isIdentical(pow(rninf, -3.0), -0.0)); 750 assert(isIdentical(pow(rninf, -2.0), +0.0)); 751 752 // @@@BUG@@@ somewhere 753 version (OSX) {} else assert(isNaN(pow(one, dinf))); 754 version (OSX) {} else assert(isNaN(pow(-one, dinf))); 755 assert(isNaN(pow(-0.2, PI))); 756 // boundary cases. Note that epsilon == 2^^-n for some n, 757 // so 1/epsilon == 2^^n is always even. 758 assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L); 759 assert(pow(-1.0L, 1/real.epsilon) == 1.0L); 760 assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L))); 761 assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L))); 762 763 assert(pow(0.0, -3.0) == double.infinity); 764 assert(pow(-0.0, -3.0) == -double.infinity); 765 assert(pow(0.0, -PI) == double.infinity); 766 assert(pow(-0.0, -PI) == double.infinity); 767 assert(isIdentical(pow(0.0, 5.0), 0.0)); 768 assert(isIdentical(pow(-0.0, 5.0), -0.0)); 769 assert(isIdentical(pow(0.0, 6.0), 0.0)); 770 assert(isIdentical(pow(-0.0, 6.0), 0.0)); 771 772 // https://issues.dlang.org/show_bug.cgi?id=14786 fixed 773 immutable real maxOdd = pow(2.0L, real.mant_dig) - 1.0L; 774 assert(pow(-1.0L, maxOdd) == -1.0L); 775 assert(pow(-1.0L, -maxOdd) == -1.0L); 776 assert(pow(-1.0L, maxOdd + 1.0L) == 1.0L); 777 assert(pow(-1.0L, -maxOdd + 1.0L) == 1.0L); 778 assert(pow(-1.0L, maxOdd - 1.0L) == 1.0L); 779 assert(pow(-1.0L, -maxOdd - 1.0L) == 1.0L); 780 781 // Now, actual numbers. 782 assert(isClose(pow(two, three), 8.0)); 783 assert(isClose(pow(two, -2.5), 0.1767766953)); 784 785 // Test integer to float power. 786 immutable uint twoI = 2; 787 assert(isClose(pow(twoI, three), 8.0)); 788 } 789 790 // https://issues.dlang.org/show_bug.cgi?id=20508 791 @safe pure nothrow @nogc unittest 792 { 793 import std.math.traits : isNaN; 794 795 assert(isNaN(pow(-double.infinity, 0.5))); 796 797 assert(isNaN(pow(-real.infinity, real.infinity))); 798 assert(isNaN(pow(-real.infinity, -real.infinity))); 799 assert(isNaN(pow(-real.infinity, 1.234))); 800 assert(isNaN(pow(-real.infinity, -0.751))); 801 assert(pow(-real.infinity, 0.0) == 1.0); 802 } 803 804 /** Computes the value of a positive integer `x`, raised to the power `n`, modulo `m`. 805 * 806 * Params: 807 * x = base 808 * n = exponent 809 * m = modulus 810 * 811 * Returns: 812 * `x` to the power `n`, modulo `m`. 813 * The return type is the largest of `x`'s and `m`'s type. 814 * 815 * The function requires that all values have unsigned types. 816 */ 817 Unqual!(Largest!(F, H)) powmod(F, G, H)(F x, G n, H m) 818 if (isUnsigned!F && isUnsigned!G && isUnsigned!H) 819 { 820 import std.meta : AliasSeq; 821 822 alias T = Unqual!(Largest!(F, H)); 823 static if (T.sizeof <= 4) 824 { 825 alias DoubleT = AliasSeq!(void, ushort, uint, void, ulong)[T.sizeof]; 826 } 827 828 static T mulmod(T a, T b, T c) 829 { 830 static if (T.sizeof == 8) 831 { 832 static T addmod(T a, T b, T c) 833 { 834 b = c - b; 835 if (a >= b) 836 return a - b; 837 else 838 return c - b + a; 839 } 840 841 T result = 0, tmp; 842 843 b %= c; 844 while (a > 0) 845 { 846 if (a & 1) 847 result = addmod(result, b, c); 848 849 a >>= 1; 850 b = addmod(b, b, c); 851 } 852 853 return result; 854 } 855 else 856 { 857 DoubleT result = cast(DoubleT) (cast(DoubleT) a * cast(DoubleT) b); 858 return result % c; 859 } 860 } 861 862 T base = x, result = 1, modulus = m; 863 Unqual!G exponent = n; 864 865 while (exponent > 0) 866 { 867 if (exponent & 1) 868 result = mulmod(result, base, modulus); 869 870 base = mulmod(base, base, modulus); 871 exponent >>= 1; 872 } 873 874 return result; 875 } 876 877 /// 878 @safe pure nothrow @nogc unittest 879 { 880 assert(powmod(1U, 10U, 3U) == 1); 881 assert(powmod(3U, 2U, 6U) == 3); 882 assert(powmod(5U, 5U, 15U) == 5); 883 assert(powmod(2U, 3U, 5U) == 3); 884 assert(powmod(2U, 4U, 5U) == 1); 885 assert(powmod(2U, 5U, 5U) == 2); 886 } 887 888 @safe pure nothrow @nogc unittest 889 { 890 ulong a = 18446744073709551615u, b = 20u, c = 18446744073709551610u; 891 assert(powmod(a, b, c) == 95367431640625u); 892 a = 100; b = 7919; c = 18446744073709551557u; 893 assert(powmod(a, b, c) == 18223853583554725198u); 894 a = 117; b = 7919; c = 18446744073709551557u; 895 assert(powmod(a, b, c) == 11493139548346411394u); 896 a = 134; b = 7919; c = 18446744073709551557u; 897 assert(powmod(a, b, c) == 10979163786734356774u); 898 a = 151; b = 7919; c = 18446744073709551557u; 899 assert(powmod(a, b, c) == 7023018419737782840u); 900 a = 168; b = 7919; c = 18446744073709551557u; 901 assert(powmod(a, b, c) == 58082701842386811u); 902 a = 185; b = 7919; c = 18446744073709551557u; 903 assert(powmod(a, b, c) == 17423478386299876798u); 904 a = 202; b = 7919; c = 18446744073709551557u; 905 assert(powmod(a, b, c) == 5522733478579799075u); 906 a = 219; b = 7919; c = 18446744073709551557u; 907 assert(powmod(a, b, c) == 15230218982491623487u); 908 a = 236; b = 7919; c = 18446744073709551557u; 909 assert(powmod(a, b, c) == 5198328724976436000u); 910 911 a = 0; b = 7919; c = 18446744073709551557u; 912 assert(powmod(a, b, c) == 0); 913 a = 123; b = 0; c = 18446744073709551557u; 914 assert(powmod(a, b, c) == 1); 915 916 immutable ulong a1 = 253, b1 = 7919, c1 = 18446744073709551557u; 917 assert(powmod(a1, b1, c1) == 3883707345459248860u); 918 919 uint x = 100 ,y = 7919, z = 1844674407u; 920 assert(powmod(x, y, z) == 1613100340u); 921 x = 134; y = 7919; z = 1844674407u; 922 assert(powmod(x, y, z) == 734956622u); 923 x = 151; y = 7919; z = 1844674407u; 924 assert(powmod(x, y, z) == 1738696945u); 925 x = 168; y = 7919; z = 1844674407u; 926 assert(powmod(x, y, z) == 1247580927u); 927 x = 185; y = 7919; z = 1844674407u; 928 assert(powmod(x, y, z) == 1293855176u); 929 x = 202; y = 7919; z = 1844674407u; 930 assert(powmod(x, y, z) == 1566963682u); 931 x = 219; y = 7919; z = 1844674407u; 932 assert(powmod(x, y, z) == 181227807u); 933 x = 236; y = 7919; z = 1844674407u; 934 assert(powmod(x, y, z) == 217988321u); 935 x = 253; y = 7919; z = 1844674407u; 936 assert(powmod(x, y, z) == 1588843243u); 937 938 x = 0; y = 7919; z = 184467u; 939 assert(powmod(x, y, z) == 0); 940 x = 123; y = 0; z = 1844674u; 941 assert(powmod(x, y, z) == 1); 942 943 immutable ubyte x1 = 117; 944 immutable uint y1 = 7919; 945 immutable uint z1 = 1844674407u; 946 auto res = powmod(x1, y1, z1); 947 assert(is(typeof(res) == uint)); 948 assert(res == 9479781u); 949 950 immutable ushort x2 = 123; 951 immutable uint y2 = 203; 952 immutable ubyte z2 = 113; 953 auto res2 = powmod(x2, y2, z2); 954 assert(is(typeof(res2) == ushort)); 955 assert(res2 == 42u); 956 } 957 958 /** 959 * Calculates e$(SUPERSCRIPT x). 960 * 961 * $(TABLE_SV 962 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)) ) 963 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) 964 * $(TR $(TD -$(INFIN)) $(TD +0.0) ) 965 * $(TR $(TD $(NAN)) $(TD $(NAN)) ) 966 * ) 967 */ 968 pragma(inline, true) 969 real exp(real x) @trusted pure nothrow @nogc // TODO: @safe 970 { 971 import std.math.constants : LOG2E; 972 973 version (InlineAsm_X87) 974 { 975 // e^^x = 2^^(LOG2E*x) 976 // (This is valid because the overflow & underflow limits for exp 977 // and exp2 are so similar). 978 if (!__ctfe) 979 return exp2Asm(LOG2E*x); 980 } 981 return expImpl(x); 982 } 983 984 /// ditto 985 pragma(inline, true) 986 double exp(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) exp(cast(real) x) : expImpl(x); } 987 988 /// ditto 989 pragma(inline, true) 990 float exp(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) exp(cast(real) x) : expImpl(x); } 991 992 /// 993 @safe unittest 994 { 995 import std.math.operations : feqrel; 996 import std.math.constants : E; 997 998 assert(exp(0.0) == 1.0); 999 assert(exp(3.0).feqrel(E * E * E) > 16); 1000 } 1001 1002 private T expImpl(T)(T x) @safe pure nothrow @nogc 1003 { 1004 import std.math : floatTraits, RealFormat; 1005 import std.math.traits : isNaN; 1006 import std.math.rounding : floor; 1007 import std.math.algebraic : poly; 1008 import std.math.constants : LOG2E; 1009 1010 alias F = floatTraits!T; 1011 static if (F.realFormat == RealFormat.ieeeSingle) 1012 { 1013 static immutable T[6] P = [ 1014 5.0000001201E-1, 1015 1.6666665459E-1, 1016 4.1665795894E-2, 1017 8.3334519073E-3, 1018 1.3981999507E-3, 1019 1.9875691500E-4, 1020 ]; 1021 1022 enum T C1 = 0.693359375; 1023 enum T C2 = -2.12194440e-4; 1024 1025 // Overflow and Underflow limits. 1026 enum T OF = 88.72283905206835; 1027 enum T UF = -103.278929903431851103; // ln(2^-149) 1028 } 1029 else static if (F.realFormat == RealFormat.ieeeDouble) 1030 { 1031 // Coefficients for exp(x) 1032 static immutable T[3] P = [ 1033 9.99999999999999999910E-1L, 1034 3.02994407707441961300E-2L, 1035 1.26177193074810590878E-4L, 1036 ]; 1037 static immutable T[4] Q = [ 1038 2.00000000000000000009E0L, 1039 2.27265548208155028766E-1L, 1040 2.52448340349684104192E-3L, 1041 3.00198505138664455042E-6L, 1042 ]; 1043 1044 // C1 + C2 = LN2. 1045 enum T C1 = 6.93145751953125E-1; 1046 enum T C2 = 1.42860682030941723212E-6; 1047 1048 // Overflow and Underflow limits. 1049 enum T OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024) 1050 enum T UF = -7.451332191019412076235E2; // ln(2^-1075) 1051 } 1052 else static if (F.realFormat == RealFormat.ieeeExtended || 1053 F.realFormat == RealFormat.ieeeExtended53) 1054 { 1055 // Coefficients for exp(x) 1056 static immutable T[3] P = [ 1057 9.9999999999999999991025E-1L, 1058 3.0299440770744196129956E-2L, 1059 1.2617719307481059087798E-4L, 1060 ]; 1061 static immutable T[4] Q = [ 1062 2.0000000000000000000897E0L, 1063 2.2726554820815502876593E-1L, 1064 2.5244834034968410419224E-3L, 1065 3.0019850513866445504159E-6L, 1066 ]; 1067 1068 // C1 + C2 = LN2. 1069 enum T C1 = 6.9314575195312500000000E-1L; 1070 enum T C2 = 1.4286068203094172321215E-6L; 1071 1072 // Overflow and Underflow limits. 1073 enum T OF = 1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384) 1074 enum T UF = -1.13994985314888605586758E4L; // ln(2^-16446) 1075 } 1076 else static if (F.realFormat == RealFormat.ieeeQuadruple) 1077 { 1078 // Coefficients for exp(x) - 1 1079 static immutable T[5] P = [ 1080 9.999999999999999999999999999999999998502E-1L, 1081 3.508710990737834361215404761139478627390E-2L, 1082 2.708775201978218837374512615596512792224E-4L, 1083 6.141506007208645008909088812338454698548E-7L, 1084 3.279723985560247033712687707263393506266E-10L 1085 ]; 1086 static immutable T[6] Q = [ 1087 2.000000000000000000000000000000000000150E0, 1088 2.368408864814233538909747618894558968880E-1L, 1089 3.611828913847589925056132680618007270344E-3L, 1090 1.504792651814944826817779302637284053660E-5L, 1091 1.771372078166251484503904874657985291164E-8L, 1092 2.980756652081995192255342779918052538681E-12L 1093 ]; 1094 1095 // C1 + C2 = LN2. 1096 enum T C1 = 6.93145751953125E-1L; 1097 enum T C2 = 1.428606820309417232121458176568075500134E-6L; 1098 1099 // Overflow and Underflow limits. 1100 enum T OF = 1.135583025911358400418251384584930671458833e4L; 1101 enum T UF = -1.143276959615573793352782661133116431383730e4L; 1102 } 1103 else 1104 static assert(0, "Not implemented for this architecture"); 1105 1106 // Special cases. 1107 if (isNaN(x)) 1108 return x; 1109 if (x > OF) 1110 return real.infinity; 1111 if (x < UF) 1112 return 0.0; 1113 1114 // Express: e^^x = e^^g * 2^^n 1115 // = e^^g * e^^(n * LOG2E) 1116 // = e^^(g + n * LOG2E) 1117 T xx = floor((cast(T) LOG2E) * x + cast(T) 0.5); 1118 const int n = cast(int) xx; 1119 x -= xx * C1; 1120 x -= xx * C2; 1121 1122 static if (F.realFormat == RealFormat.ieeeSingle) 1123 { 1124 xx = x * x; 1125 x = poly(x, P) * xx + x + 1.0f; 1126 } 1127 else 1128 { 1129 // Rational approximation for exponential of the fractional part: 1130 // e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2)) 1131 xx = x * x; 1132 const T px = x * poly(xx, P); 1133 x = px / (poly(xx, Q) - px); 1134 x = (cast(T) 1.0) + (cast(T) 2.0) * x; 1135 } 1136 1137 // Scale by power of 2. 1138 x = core.math.ldexp(x, n); 1139 1140 return x; 1141 } 1142 1143 @safe @nogc nothrow unittest 1144 { 1145 import std.math : floatTraits, RealFormat; 1146 import std.math.operations : NaN, feqrel, isClose; 1147 import std.math.constants : E; 1148 import std.math.traits : isIdentical; 1149 import std.math.algebraic : abs; 1150 1151 version (IeeeFlagsSupport) import std.math.hardware : IeeeFlags, resetIeeeFlags, ieeeFlags; 1152 version (FloatingPointControlSupport) 1153 { 1154 import std.math.hardware : FloatingPointControl; 1155 1156 FloatingPointControl ctrl; 1157 if (FloatingPointControl.hasExceptionTraps) 1158 ctrl.disableExceptions(FloatingPointControl.allExceptions); 1159 ctrl.rounding = FloatingPointControl.roundToNearest; 1160 } 1161 1162 static void testExp(T)() 1163 { 1164 enum realFormat = floatTraits!T.realFormat; 1165 static if (realFormat == RealFormat.ieeeQuadruple) 1166 { 1167 static immutable T[2][] exptestpoints = 1168 [ // x exp(x) 1169 [ 1.0L, E ], 1170 [ 0.5L, 0x1.a61298e1e069bc972dfefab6df34p+0L ], 1171 [ 3.0L, E*E*E ], 1172 [ 0x1.6p+13L, 0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow 1173 [ 0x1.7p+13L, T.infinity ], // close overflow 1174 [ 0x1p+80L, T.infinity ], // far overflow 1175 [ T.infinity, T.infinity ], 1176 [-0x1.18p+13L, 0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow 1177 [-0x1.625p+13L, 0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto 1178 [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal 1179 [-0x1.6549p+13L, 0x0.0000000000000000000000000001p-16382L ], // ditto 1180 [-0x1.655p+13L, 0 ], // close underflow 1181 [-0x1p+30L, 0 ], // far underflow 1182 ]; 1183 } 1184 else static if (realFormat == RealFormat.ieeeExtended || 1185 realFormat == RealFormat.ieeeExtended53) 1186 { 1187 static immutable T[2][] exptestpoints = 1188 [ // x exp(x) 1189 [ 1.0L, E ], 1190 [ 0.5L, 0x1.a61298e1e069bc97p+0L ], 1191 [ 3.0L, E*E*E ], 1192 [ 0x1.1p+13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow 1193 [ 0x1.7p+13L, T.infinity ], // close overflow 1194 [ 0x1p+80L, T.infinity ], // far overflow 1195 [ T.infinity, T.infinity ], 1196 [-0x1.18p+13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow 1197 [-0x1.625p+13L, 0x1.a6bd68a39d11f35cp-16358L ], // ditto 1198 [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L ], // near underflow - subnormal 1199 [-0x1.643p+13L, 0x1p-16444L ], // ditto 1200 [-0x1.645p+13L, 0 ], // close underflow 1201 [-0x1p+30L, 0 ], // far underflow 1202 ]; 1203 } 1204 else static if (realFormat == RealFormat.ieeeDouble) 1205 { 1206 static immutable T[2][] exptestpoints = 1207 [ // x, exp(x) 1208 [ 1.0L, E ], 1209 [ 0.5L, 0x1.a61298e1e069cp+0L ], 1210 [ 3.0L, E*E*E ], 1211 [ 0x1.6p+9L, 0x1.93bf4ec282efbp+1015L ], // near overflow 1212 [ 0x1.7p+9L, T.infinity ], // close overflow 1213 [ 0x1p+80L, T.infinity ], // far overflow 1214 [ T.infinity, T.infinity ], 1215 [-0x1.6p+9L, 0x1.44a3824e5285fp-1016L ], // near underflow 1216 [-0x1.64p+9L, 0x0.06f84920bb2d4p-1022L ], // near underflow - subnormal 1217 [-0x1.743p+9L, 0x0.0000000000001p-1022L ], // ditto 1218 [-0x1.8p+9L, 0 ], // close underflow 1219 [-0x1p+30L, 0 ], // far underflow 1220 ]; 1221 } 1222 else static if (realFormat == RealFormat.ieeeSingle) 1223 { 1224 static immutable T[2][] exptestpoints = 1225 [ // x, exp(x) 1226 [ 1.0L, E ], 1227 [ 0.5L, 0x1.a61299p+0L ], 1228 [ 3.0L, E*E*E ], 1229 [ 0x1.62p+6L, 0x1.99b988p+127L ], // near overflow 1230 [ 0x1.7p+6L, T.infinity ], // close overflow 1231 [ 0x1p+80L, T.infinity ], // far overflow 1232 [ T.infinity, T.infinity ], 1233 [-0x1.5cp+6L, 0x1.666d0ep-126L ], // near underflow 1234 [-0x1.7p+6L, 0x0.026a42p-126L ], // near underflow - subnormal 1235 [-0x1.9cp+6L, 0x0.000002p-126L ], // ditto 1236 [-0x1.ap+6L, 0 ], // close underflow 1237 [-0x1p+30L, 0 ], // far underflow 1238 ]; 1239 } 1240 else 1241 static assert(0, "No exp() tests for real type!"); 1242 1243 const minEqualMantissaBits = T.mant_dig - 13; 1244 T x; 1245 version (IeeeFlagsSupport) IeeeFlags f; 1246 foreach (ref pair; exptestpoints) 1247 { 1248 version (IeeeFlagsSupport) resetIeeeFlags(); 1249 x = exp(pair[0]); 1250 //printf("exp(%La) = %La, should be %La\n", cast(real) pair[0], cast(real) x, cast(real) pair[1]); 1251 assert(feqrel(x, pair[1]) >= minEqualMantissaBits); 1252 } 1253 1254 // Ideally, exp(0) would not set the inexact flag. 1255 // Unfortunately, fldl2e sets it! 1256 // So it's not realistic to avoid setting it. 1257 assert(exp(cast(T) 0.0) == 1.0); 1258 1259 // NaN propagation. Doesn't set flags, bcos was already NaN. 1260 version (IeeeFlagsSupport) 1261 { 1262 resetIeeeFlags(); 1263 x = exp(T.nan); 1264 f = ieeeFlags; 1265 assert(isIdentical(abs(x), T.nan)); 1266 assert(f.flags == 0); 1267 1268 resetIeeeFlags(); 1269 x = exp(-T.nan); 1270 f = ieeeFlags; 1271 assert(isIdentical(abs(x), T.nan)); 1272 assert(f.flags == 0); 1273 } 1274 else 1275 { 1276 x = exp(T.nan); 1277 assert(isIdentical(abs(x), T.nan)); 1278 1279 x = exp(-T.nan); 1280 assert(isIdentical(abs(x), T.nan)); 1281 } 1282 1283 x = exp(NaN(0x123)); 1284 assert(isIdentical(x, NaN(0x123))); 1285 } 1286 1287 import std.meta : AliasSeq; 1288 foreach (T; AliasSeq!(real, double, float)) 1289 testExp!T(); 1290 1291 // High resolution test (verified against GNU MPFR/Mathematica). 1292 assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+0L); 1293 1294 assert(isClose(exp(3.0L), E * E * E, real.sizeof > double.sizeof ? 1e-15 : 1e-14)); 1295 } 1296 1297 /** 1298 * Calculates the value of the natural logarithm base (e) 1299 * raised to the power of x, minus 1. 1300 * 1301 * For very small x, expm1(x) is more accurate 1302 * than exp(x)-1. 1303 * 1304 * $(TABLE_SV 1305 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)-1) ) 1306 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) ) 1307 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) 1308 * $(TR $(TD -$(INFIN)) $(TD -1.0) ) 1309 * $(TR $(TD $(NAN)) $(TD $(NAN)) ) 1310 * ) 1311 */ 1312 pragma(inline, true) 1313 real expm1(real x) @trusted pure nothrow @nogc // TODO: @safe 1314 { 1315 version (InlineAsm_X87) 1316 { 1317 if (!__ctfe) 1318 return expm1Asm(x); 1319 } 1320 return expm1Impl(x); 1321 } 1322 1323 /// ditto 1324 pragma(inline, true) 1325 double expm1(double x) @safe pure nothrow @nogc 1326 { 1327 return __ctfe ? cast(double) expm1(cast(real) x) : expm1Impl(x); 1328 } 1329 1330 /// ditto 1331 pragma(inline, true) 1332 float expm1(float x) @safe pure nothrow @nogc 1333 { 1334 // no single-precision version in Cephes => use double precision 1335 return __ctfe ? cast(float) expm1(cast(real) x) : cast(float) expm1Impl(cast(double) x); 1336 } 1337 1338 /// 1339 @safe unittest 1340 { 1341 import std.math.traits : isIdentical; 1342 import std.math.operations : feqrel; 1343 1344 assert(isIdentical(expm1(0.0), 0.0)); 1345 assert(expm1(1.0).feqrel(1.71828) > 16); 1346 assert(expm1(2.0).feqrel(6.3890) > 16); 1347 } 1348 1349 version (InlineAsm_X87) 1350 private real expm1Asm(real x) @trusted pure nothrow @nogc 1351 { 1352 version (X86) 1353 { 1354 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4 1355 asm pure nothrow @nogc 1356 { 1357 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. 1358 * Author: Don Clugston. 1359 * 1360 * expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x. 1361 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y)) 1362 * and 2ym1 = (2^^(y-rndint(y))-1). 1363 * If 2rndy < 0.5*real.epsilon, result is -1. 1364 * Implementation is otherwise the same as for exp2() 1365 */ 1366 naked; 1367 fld real ptr [ESP+4] ; // x 1368 mov AX, [ESP+4+8]; // AX = exponent and sign 1369 sub ESP, 12+8; // Create scratch space on the stack 1370 // [ESP,ESP+2] = scratchint 1371 // [ESP+4..+6, +8..+10, +10] = scratchreal 1372 // set scratchreal mantissa = 1.0 1373 mov dword ptr [ESP+8], 0; 1374 mov dword ptr [ESP+8+4], 0x80000000; 1375 and AX, 0x7FFF; // drop sign bit 1376 cmp AX, 0x401D; // avoid InvalidException in fist 1377 jae L_extreme; 1378 fldl2e; 1379 fmulp ST(1), ST; // y = x*log2(e) 1380 fist dword ptr [ESP]; // scratchint = rndint(y) 1381 fisub dword ptr [ESP]; // y - rndint(y) 1382 // and now set scratchreal exponent 1383 mov EAX, [ESP]; 1384 add EAX, 0x3fff; 1385 jle short L_largenegative; 1386 cmp EAX,0x8000; 1387 jge short L_largepositive; 1388 mov [ESP+8+8],AX; 1389 f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1 1390 fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y) 1391 fmul ST(1), ST; // ST=2rndy, ST(1)=2rndy*2ym1 1392 fld1; 1393 fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1 1394 faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1 1395 add ESP,12+8; 1396 ret PARAMSIZE; 1397 1398 L_extreme: // Extreme exponent. X is very large positive, very 1399 // large negative, infinity, or NaN. 1400 fxam; 1401 fstsw AX; 1402 test AX, 0x0400; // NaN_or_zero, but we already know x != 0 1403 jz L_was_nan; // if x is NaN, returns x 1404 test AX, 0x0200; 1405 jnz L_largenegative; 1406 L_largepositive: 1407 // Set scratchreal = real.max. 1408 // squaring it will create infinity, and set overflow flag. 1409 mov word ptr [ESP+8+8], 0x7FFE; 1410 fstp ST(0); 1411 fld real ptr [ESP+8]; // load scratchreal 1412 fmul ST(0), ST; // square it, to create havoc! 1413 L_was_nan: 1414 add ESP,12+8; 1415 ret PARAMSIZE; 1416 L_largenegative: 1417 fstp ST(0); 1418 fld1; 1419 fchs; // return -1. Underflow flag is not set. 1420 add ESP,12+8; 1421 ret PARAMSIZE; 1422 } 1423 } 1424 else version (X86_64) 1425 { 1426 asm pure nothrow @nogc 1427 { 1428 naked; 1429 } 1430 version (Win64) 1431 { 1432 asm pure nothrow @nogc 1433 { 1434 fld real ptr [RCX]; // x 1435 mov AX,[RCX+8]; // AX = exponent and sign 1436 } 1437 } 1438 else 1439 { 1440 asm pure nothrow @nogc 1441 { 1442 fld real ptr [RSP+8]; // x 1443 mov AX,[RSP+8+8]; // AX = exponent and sign 1444 } 1445 } 1446 asm pure nothrow @nogc 1447 { 1448 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. 1449 * Author: Don Clugston. 1450 * 1451 * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x. 1452 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y)) 1453 * and 2ym1 = (2^(y-rndint(y))-1). 1454 * If 2rndy < 0.5*real.epsilon, result is -1. 1455 * Implementation is otherwise the same as for exp2() 1456 */ 1457 sub RSP, 24; // Create scratch space on the stack 1458 // [RSP,RSP+2] = scratchint 1459 // [RSP+4..+6, +8..+10, +10] = scratchreal 1460 // set scratchreal mantissa = 1.0 1461 mov dword ptr [RSP+8], 0; 1462 mov dword ptr [RSP+8+4], 0x80000000; 1463 and AX, 0x7FFF; // drop sign bit 1464 cmp AX, 0x401D; // avoid InvalidException in fist 1465 jae L_extreme; 1466 fldl2e; 1467 fmul ; // y = x*log2(e) 1468 fist dword ptr [RSP]; // scratchint = rndint(y) 1469 fisub dword ptr [RSP]; // y - rndint(y) 1470 // and now set scratchreal exponent 1471 mov EAX, [RSP]; 1472 add EAX, 0x3fff; 1473 jle short L_largenegative; 1474 cmp EAX,0x8000; 1475 jge short L_largepositive; 1476 mov [RSP+8+8],AX; 1477 f2xm1; // 2^(y-rndint(y)) -1 1478 fld real ptr [RSP+8] ; // 2^rndint(y) 1479 fmul ST(1), ST; 1480 fld1; 1481 fsubp ST(1), ST; 1482 fadd; 1483 add RSP,24; 1484 ret; 1485 1486 L_extreme: // Extreme exponent. X is very large positive, very 1487 // large negative, infinity, or NaN. 1488 fxam; 1489 fstsw AX; 1490 test AX, 0x0400; // NaN_or_zero, but we already know x != 0 1491 jz L_was_nan; // if x is NaN, returns x 1492 test AX, 0x0200; 1493 jnz L_largenegative; 1494 L_largepositive: 1495 // Set scratchreal = real.max. 1496 // squaring it will create infinity, and set overflow flag. 1497 mov word ptr [RSP+8+8], 0x7FFE; 1498 fstp ST(0); 1499 fld real ptr [RSP+8]; // load scratchreal 1500 fmul ST(0), ST; // square it, to create havoc! 1501 L_was_nan: 1502 add RSP,24; 1503 ret; 1504 1505 L_largenegative: 1506 fstp ST(0); 1507 fld1; 1508 fchs; // return -1. Underflow flag is not set. 1509 add RSP,24; 1510 ret; 1511 } 1512 } 1513 else 1514 static assert(0); 1515 } 1516 1517 private T expm1Impl(T)(T x) @safe pure nothrow @nogc 1518 { 1519 import std.math : floatTraits, RealFormat; 1520 import std.math.rounding : floor; 1521 import std.math.algebraic : poly; 1522 import std.math.constants : LN2; 1523 1524 // Coefficients for exp(x) - 1 and overflow/underflow limits. 1525 enum realFormat = floatTraits!T.realFormat; 1526 static if (realFormat == RealFormat.ieeeQuadruple) 1527 { 1528 static immutable T[8] P = [ 1529 2.943520915569954073888921213330863757240E8L, 1530 -5.722847283900608941516165725053359168840E7L, 1531 8.944630806357575461578107295909719817253E6L, 1532 -7.212432713558031519943281748462837065308E5L, 1533 4.578962475841642634225390068461943438441E4L, 1534 -1.716772506388927649032068540558788106762E3L, 1535 4.401308817383362136048032038528753151144E1L, 1536 -4.888737542888633647784737721812546636240E-1L 1537 ]; 1538 1539 static immutable T[9] Q = [ 1540 1.766112549341972444333352727998584753865E9L, 1541 -7.848989743695296475743081255027098295771E8L, 1542 1.615869009634292424463780387327037251069E8L, 1543 -2.019684072836541751428967854947019415698E7L, 1544 1.682912729190313538934190635536631941751E6L, 1545 -9.615511549171441430850103489315371768998E4L, 1546 3.697714952261803935521187272204485251835E3L, 1547 -8.802340681794263968892934703309274564037E1L, 1548 1.0 1549 ]; 1550 1551 enum T OF = 1.1356523406294143949491931077970764891253E4L; 1552 enum T UF = -1.143276959615573793352782661133116431383730e4L; 1553 } 1554 else static if (realFormat == RealFormat.ieeeExtended) 1555 { 1556 static immutable T[5] P = [ 1557 -1.586135578666346600772998894928250240826E4L, 1558 2.642771505685952966904660652518429479531E3L, 1559 -3.423199068835684263987132888286791620673E2L, 1560 1.800826371455042224581246202420972737840E1L, 1561 -5.238523121205561042771939008061958820811E-1L, 1562 ]; 1563 static immutable T[6] Q = [ 1564 -9.516813471998079611319047060563358064497E4L, 1565 3.964866271411091674556850458227710004570E4L, 1566 -7.207678383830091850230366618190187434796E3L, 1567 7.206038318724600171970199625081491823079E2L, 1568 -4.002027679107076077238836622982900945173E1L, 1569 1.0 1570 ]; 1571 1572 enum T OF = 1.1356523406294143949492E4L; 1573 enum T UF = -4.5054566736396445112120088E1L; 1574 } 1575 else static if (realFormat == RealFormat.ieeeDouble) 1576 { 1577 static immutable T[3] P = [ 1578 9.9999999999999999991025E-1, 1579 3.0299440770744196129956E-2, 1580 1.2617719307481059087798E-4, 1581 ]; 1582 static immutable T[4] Q = [ 1583 2.0000000000000000000897E0, 1584 2.2726554820815502876593E-1, 1585 2.5244834034968410419224E-3, 1586 3.0019850513866445504159E-6, 1587 ]; 1588 } 1589 else 1590 static assert(0, "no coefficients for expm1()"); 1591 1592 static if (realFormat == RealFormat.ieeeDouble) // special case for double precision 1593 { 1594 if (x < -0.5 || x > 0.5) 1595 return exp(x) - 1.0; 1596 if (x == 0.0) 1597 return x; 1598 1599 const T xx = x * x; 1600 x = x * poly(xx, P); 1601 x = x / (poly(xx, Q) - x); 1602 return x + x; 1603 } 1604 else 1605 { 1606 // C1 + C2 = LN2. 1607 enum T C1 = 6.9314575195312500000000E-1L; 1608 enum T C2 = 1.428606820309417232121458176568075500134E-6L; 1609 1610 // Special cases. 1611 if (x > OF) 1612 return real.infinity; 1613 if (x == cast(T) 0.0) 1614 return x; 1615 if (x < UF) 1616 return -1.0; 1617 1618 // Express x = LN2 (n + remainder), remainder not exceeding 1/2. 1619 int n = cast(int) floor((cast(T) 0.5) + x / cast(T) LN2); 1620 x -= n * C1; 1621 x -= n * C2; 1622 1623 // Rational approximation: 1624 // exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x) 1625 T px = x * poly(x, P); 1626 T qx = poly(x, Q); 1627 const T xx = x * x; 1628 qx = x + ((cast(T) 0.5) * xx + xx * px / qx); 1629 1630 // We have qx = exp(remainder LN2) - 1, so: 1631 // exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1. 1632 px = core.math.ldexp(cast(T) 1.0, n); 1633 x = px * qx + (px - cast(T) 1.0); 1634 1635 return x; 1636 } 1637 } 1638 1639 @safe @nogc nothrow unittest 1640 { 1641 import std.math.traits : isNaN; 1642 import std.math.operations : isClose, CommonDefaultFor; 1643 1644 static void testExpm1(T)() 1645 { 1646 // NaN 1647 assert(isNaN(expm1(cast(T) T.nan))); 1648 1649 static immutable T[] xs = [ -2, -0.75, -0.3, 0.0, 0.1, 0.2, 0.5, 1.0 ]; 1650 foreach (x; xs) 1651 { 1652 const T e = expm1(x); 1653 const T r = exp(x) - 1; 1654 1655 //printf("expm1(%Lg) = %Lg, should approximately be %Lg\n", cast(real) x, cast(real) e, cast(real) r); 1656 assert(isClose(r, e, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T))); 1657 } 1658 } 1659 1660 import std.meta : AliasSeq; 1661 foreach (T; AliasSeq!(real, double)) 1662 testExpm1!T(); 1663 } 1664 1665 /** 1666 * Calculates 2$(SUPERSCRIPT x). 1667 * 1668 * $(TABLE_SV 1669 * $(TR $(TH x) $(TH exp2(x)) ) 1670 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) 1671 * $(TR $(TD -$(INFIN)) $(TD +0.0) ) 1672 * $(TR $(TD $(NAN)) $(TD $(NAN)) ) 1673 * ) 1674 */ 1675 pragma(inline, true) 1676 real exp2(real x) @nogc @trusted pure nothrow // TODO: @safe 1677 { 1678 version (InlineAsm_X87) 1679 { 1680 if (!__ctfe) 1681 return exp2Asm(x); 1682 } 1683 return exp2Impl(x); 1684 } 1685 1686 /// ditto 1687 pragma(inline, true) 1688 double exp2(double x) @nogc @safe pure nothrow { return __ctfe ? cast(double) exp2(cast(real) x) : exp2Impl(x); } 1689 1690 /// ditto 1691 pragma(inline, true) 1692 float exp2(float x) @nogc @safe pure nothrow { return __ctfe ? cast(float) exp2(cast(real) x) : exp2Impl(x); } 1693 1694 /// 1695 @safe unittest 1696 { 1697 import std.math.traits : isIdentical; 1698 import std.math.operations : feqrel; 1699 1700 assert(isIdentical(exp2(0.0), 1.0)); 1701 assert(exp2(2.0).feqrel(4.0) > 16); 1702 assert(exp2(8.0).feqrel(256.0) > 16); 1703 } 1704 1705 @safe unittest 1706 { 1707 version (CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented 1708 { 1709 assert( core.stdc.math.exp2f(0.0f) == 1 ); 1710 assert( core.stdc.math.exp2 (0.0) == 1 ); 1711 assert( core.stdc.math.exp2l(0.0L) == 1 ); 1712 } 1713 } 1714 1715 version (InlineAsm_X87) 1716 private real exp2Asm(real x) @nogc @trusted pure nothrow 1717 { 1718 version (X86) 1719 { 1720 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4 1721 1722 asm pure nothrow @nogc 1723 { 1724 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. 1725 * Author: Don Clugston. 1726 * 1727 * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x)) 1728 * The trick for high performance is to avoid the fscale(28cycles on core2), 1729 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. 1730 * 1731 * We can do frndint by using fist. BUT we can't use it for huge numbers, 1732 * because it will set the Invalid Operation flag if overflow or NaN occurs. 1733 * Fortunately, whenever this happens the result would be zero or infinity. 1734 * 1735 * We can perform fscale by directly poking into the exponent. BUT this doesn't 1736 * work for the (very rare) cases where the result is subnormal. So we fall back 1737 * to the slow method in that case. 1738 */ 1739 naked; 1740 fld real ptr [ESP+4] ; // x 1741 mov AX, [ESP+4+8]; // AX = exponent and sign 1742 sub ESP, 12+8; // Create scratch space on the stack 1743 // [ESP,ESP+2] = scratchint 1744 // [ESP+4..+6, +8..+10, +10] = scratchreal 1745 // set scratchreal mantissa = 1.0 1746 mov dword ptr [ESP+8], 0; 1747 mov dword ptr [ESP+8+4], 0x80000000; 1748 and AX, 0x7FFF; // drop sign bit 1749 cmp AX, 0x401D; // avoid InvalidException in fist 1750 jae L_extreme; 1751 fist dword ptr [ESP]; // scratchint = rndint(x) 1752 fisub dword ptr [ESP]; // x - rndint(x) 1753 // and now set scratchreal exponent 1754 mov EAX, [ESP]; 1755 add EAX, 0x3fff; 1756 jle short L_subnormal; 1757 cmp EAX,0x8000; 1758 jge short L_overflow; 1759 mov [ESP+8+8],AX; 1760 L_normal: 1761 f2xm1; 1762 fld1; 1763 faddp ST(1), ST; // 2^^(x-rndint(x)) 1764 fld real ptr [ESP+8] ; // 2^^rndint(x) 1765 add ESP,12+8; 1766 fmulp ST(1), ST; 1767 ret PARAMSIZE; 1768 1769 L_subnormal: 1770 // Result will be subnormal. 1771 // In this rare case, the simple poking method doesn't work. 1772 // The speed doesn't matter, so use the slow fscale method. 1773 fild dword ptr [ESP]; // scratchint 1774 fld1; 1775 fscale; 1776 fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint 1777 fstp ST(0); // drop scratchint 1778 jmp L_normal; 1779 1780 L_extreme: // Extreme exponent. X is very large positive, very 1781 // large negative, infinity, or NaN. 1782 fxam; 1783 fstsw AX; 1784 test AX, 0x0400; // NaN_or_zero, but we already know x != 0 1785 jz L_was_nan; // if x is NaN, returns x 1786 // set scratchreal = real.min_normal 1787 // squaring it will return 0, setting underflow flag 1788 mov word ptr [ESP+8+8], 1; 1789 test AX, 0x0200; 1790 jnz L_waslargenegative; 1791 L_overflow: 1792 // Set scratchreal = real.max. 1793 // squaring it will create infinity, and set overflow flag. 1794 mov word ptr [ESP+8+8], 0x7FFE; 1795 L_waslargenegative: 1796 fstp ST(0); 1797 fld real ptr [ESP+8]; // load scratchreal 1798 fmul ST(0), ST; // square it, to create havoc! 1799 L_was_nan: 1800 add ESP,12+8; 1801 ret PARAMSIZE; 1802 } 1803 } 1804 else version (X86_64) 1805 { 1806 asm pure nothrow @nogc 1807 { 1808 naked; 1809 } 1810 version (Win64) 1811 { 1812 asm pure nothrow @nogc 1813 { 1814 fld real ptr [RCX]; // x 1815 mov AX,[RCX+8]; // AX = exponent and sign 1816 } 1817 } 1818 else 1819 { 1820 asm pure nothrow @nogc 1821 { 1822 fld real ptr [RSP+8]; // x 1823 mov AX,[RSP+8+8]; // AX = exponent and sign 1824 } 1825 } 1826 asm pure nothrow @nogc 1827 { 1828 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. 1829 * Author: Don Clugston. 1830 * 1831 * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x)) 1832 * The trick for high performance is to avoid the fscale(28cycles on core2), 1833 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. 1834 * 1835 * We can do frndint by using fist. BUT we can't use it for huge numbers, 1836 * because it will set the Invalid Operation flag is overflow or NaN occurs. 1837 * Fortunately, whenever this happens the result would be zero or infinity. 1838 * 1839 * We can perform fscale by directly poking into the exponent. BUT this doesn't 1840 * work for the (very rare) cases where the result is subnormal. So we fall back 1841 * to the slow method in that case. 1842 */ 1843 sub RSP, 24; // Create scratch space on the stack 1844 // [RSP,RSP+2] = scratchint 1845 // [RSP+4..+6, +8..+10, +10] = scratchreal 1846 // set scratchreal mantissa = 1.0 1847 mov dword ptr [RSP+8], 0; 1848 mov dword ptr [RSP+8+4], 0x80000000; 1849 and AX, 0x7FFF; // drop sign bit 1850 cmp AX, 0x401D; // avoid InvalidException in fist 1851 jae L_extreme; 1852 fist dword ptr [RSP]; // scratchint = rndint(x) 1853 fisub dword ptr [RSP]; // x - rndint(x) 1854 // and now set scratchreal exponent 1855 mov EAX, [RSP]; 1856 add EAX, 0x3fff; 1857 jle short L_subnormal; 1858 cmp EAX,0x8000; 1859 jge short L_overflow; 1860 mov [RSP+8+8],AX; 1861 L_normal: 1862 f2xm1; 1863 fld1; 1864 fadd; // 2^(x-rndint(x)) 1865 fld real ptr [RSP+8] ; // 2^rndint(x) 1866 add RSP,24; 1867 fmulp ST(1), ST; 1868 ret; 1869 1870 L_subnormal: 1871 // Result will be subnormal. 1872 // In this rare case, the simple poking method doesn't work. 1873 // The speed doesn't matter, so use the slow fscale method. 1874 fild dword ptr [RSP]; // scratchint 1875 fld1; 1876 fscale; 1877 fstp real ptr [RSP+8]; // scratchreal = 2^scratchint 1878 fstp ST(0); // drop scratchint 1879 jmp L_normal; 1880 1881 L_extreme: // Extreme exponent. X is very large positive, very 1882 // large negative, infinity, or NaN. 1883 fxam; 1884 fstsw AX; 1885 test AX, 0x0400; // NaN_or_zero, but we already know x != 0 1886 jz L_was_nan; // if x is NaN, returns x 1887 // set scratchreal = real.min 1888 // squaring it will return 0, setting underflow flag 1889 mov word ptr [RSP+8+8], 1; 1890 test AX, 0x0200; 1891 jnz L_waslargenegative; 1892 L_overflow: 1893 // Set scratchreal = real.max. 1894 // squaring it will create infinity, and set overflow flag. 1895 mov word ptr [RSP+8+8], 0x7FFE; 1896 L_waslargenegative: 1897 fstp ST(0); 1898 fld real ptr [RSP+8]; // load scratchreal 1899 fmul ST(0), ST; // square it, to create havoc! 1900 L_was_nan: 1901 add RSP,24; 1902 ret; 1903 } 1904 } 1905 else 1906 static assert(0); 1907 } 1908 1909 private T exp2Impl(T)(T x) @nogc @safe pure nothrow 1910 { 1911 import std.math : floatTraits, RealFormat; 1912 import std.math.traits : isNaN; 1913 import std.math.rounding : floor; 1914 import std.math.algebraic : poly; 1915 1916 // Coefficients for exp2(x) 1917 enum realFormat = floatTraits!T.realFormat; 1918 static if (realFormat == RealFormat.ieeeQuadruple) 1919 { 1920 static immutable T[5] P = [ 1921 9.079594442980146270952372234833529694788E12L, 1922 1.530625323728429161131811299626419117557E11L, 1923 5.677513871931844661829755443994214173883E8L, 1924 6.185032670011643762127954396427045467506E5L, 1925 1.587171580015525194694938306936721666031E2L 1926 ]; 1927 1928 static immutable T[6] Q = [ 1929 2.619817175234089411411070339065679229869E13L, 1930 1.490560994263653042761789432690793026977E12L, 1931 1.092141473886177435056423606755843616331E10L, 1932 2.186249607051644894762167991800811827835E7L, 1933 1.236602014442099053716561665053645270207E4L, 1934 1.0 1935 ]; 1936 } 1937 else static if (realFormat == RealFormat.ieeeExtended) 1938 { 1939 static immutable T[3] P = [ 1940 2.0803843631901852422887E6L, 1941 3.0286971917562792508623E4L, 1942 6.0614853552242266094567E1L, 1943 ]; 1944 static immutable T[4] Q = [ 1945 6.0027204078348487957118E6L, 1946 3.2772515434906797273099E5L, 1947 1.7492876999891839021063E3L, 1948 1.0000000000000000000000E0L, 1949 ]; 1950 } 1951 else static if (realFormat == RealFormat.ieeeDouble) 1952 { 1953 static immutable T[3] P = [ 1954 1.51390680115615096133E3L, 1955 2.02020656693165307700E1L, 1956 2.30933477057345225087E-2L, 1957 ]; 1958 static immutable T[3] Q = [ 1959 4.36821166879210612817E3L, 1960 2.33184211722314911771E2L, 1961 1.00000000000000000000E0L, 1962 ]; 1963 } 1964 else static if (realFormat == RealFormat.ieeeSingle) 1965 { 1966 static immutable T[6] P = [ 1967 6.931472028550421E-001L, 1968 2.402264791363012E-001L, 1969 5.550332471162809E-002L, 1970 9.618437357674640E-003L, 1971 1.339887440266574E-003L, 1972 1.535336188319500E-004L, 1973 ]; 1974 } 1975 else 1976 static assert(0, "no coefficients for exp2()"); 1977 1978 // Overflow and Underflow limits. 1979 enum T OF = T.max_exp; 1980 enum T UF = T.min_exp - 1; 1981 1982 // Special cases. 1983 if (isNaN(x)) 1984 return x; 1985 if (x > OF) 1986 return real.infinity; 1987 if (x < UF) 1988 return 0.0; 1989 1990 static if (realFormat == RealFormat.ieeeSingle) // special case for single precision 1991 { 1992 // The following is necessary because range reduction blows up. 1993 if (x == 0.0f) 1994 return 1.0f; 1995 1996 // Separate into integer and fractional parts. 1997 const T i = floor(x); 1998 int n = cast(int) i; 1999 x -= i; 2000 if (x > 0.5f) 2001 { 2002 n += 1; 2003 x -= 1.0f; 2004 } 2005 2006 // Rational approximation: 2007 // exp2(x) = 1.0 + x P(x) 2008 x = 1.0f + x * poly(x, P); 2009 } 2010 else 2011 { 2012 // Separate into integer and fractional parts. 2013 const T i = floor(x + cast(T) 0.5); 2014 int n = cast(int) i; 2015 x -= i; 2016 2017 // Rational approximation: 2018 // exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2)) 2019 const T xx = x * x; 2020 const T px = x * poly(xx, P); 2021 x = px / (poly(xx, Q) - px); 2022 x = (cast(T) 1.0) + (cast(T) 2.0) * x; 2023 } 2024 2025 // Scale by power of 2. 2026 x = core.math.ldexp(x, n); 2027 2028 return x; 2029 } 2030 2031 @safe @nogc nothrow unittest 2032 { 2033 import std.math.operations : feqrel, NaN, isClose; 2034 import std.math.traits : isIdentical; 2035 import std.math.constants : SQRT2; 2036 2037 assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1); 2038 assert(exp2(8.0L) == 256.0); 2039 assert(exp2(-9.0L)== 1.0L/512.0); 2040 2041 static void testExp2(T)() 2042 { 2043 // NaN 2044 const T specialNaN = NaN(0x0123L); 2045 assert(isIdentical(exp2(specialNaN), specialNaN)); 2046 2047 // over-/underflow 2048 enum T OF = T.max_exp; 2049 enum T UF = T.min_exp - T.mant_dig; 2050 assert(isIdentical(exp2(OF + 1), cast(T) T.infinity)); 2051 assert(isIdentical(exp2(UF - 1), cast(T) 0.0)); 2052 2053 static immutable T[2][] vals = 2054 [ 2055 // x, exp2(x) 2056 [ 0.0, 1.0 ], 2057 [ -0.0, 1.0 ], 2058 [ 0.5, SQRT2 ], 2059 [ 8.0, 256.0 ], 2060 [ -9.0, 1.0 / 512 ], 2061 ]; 2062 2063 foreach (ref val; vals) 2064 { 2065 const T x = val[0]; 2066 const T r = val[1]; 2067 const T e = exp2(x); 2068 2069 //printf("exp2(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) e, cast(real) r); 2070 assert(isClose(r, e)); 2071 } 2072 } 2073 2074 import std.meta : AliasSeq; 2075 foreach (T; AliasSeq!(real, double, float)) 2076 testExp2!T(); 2077 } 2078 2079 /********************************************************************* 2080 * Separate floating point value into significand and exponent. 2081 * 2082 * Returns: 2083 * Calculate and return $(I x) and $(I exp) such that 2084 * value =$(I x)*2$(SUPERSCRIPT exp) and 2085 * .5 $(LT)= |$(I x)| $(LT) 1.0 2086 * 2087 * $(I x) has same sign as value. 2088 * 2089 * $(TABLE_SV 2090 * $(TR $(TH value) $(TH returns) $(TH exp)) 2091 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0)) 2092 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max)) 2093 * $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min)) 2094 * $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min)) 2095 * ) 2096 */ 2097 T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc 2098 if (isFloatingPoint!T) 2099 { 2100 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB; 2101 import std.math.traits : isSubnormal; 2102 2103 if (__ctfe) 2104 { 2105 // Handle special cases. 2106 if (value == 0) { exp = 0; return value; } 2107 if (value == T.infinity) { exp = int.max; return value; } 2108 if (value == -T.infinity || value != value) { exp = int.min; return value; } 2109 // Handle ordinary cases. 2110 // In CTFE there is no performance advantage for having separate 2111 // paths for different floating point types. 2112 T absValue = value < 0 ? -value : value; 2113 int expCount; 2114 static if (T.mant_dig > double.mant_dig) 2115 { 2116 for (; absValue >= 0x1.0p+1024L; absValue *= 0x1.0p-1024L) 2117 expCount += 1024; 2118 for (; absValue < 0x1.0p-1021L; absValue *= 0x1.0p+1021L) 2119 expCount -= 1021; 2120 } 2121 const double dval = cast(double) absValue; 2122 int dexp = cast(int) (((*cast(const long*) &dval) >>> 52) & 0x7FF) + double.min_exp - 2; 2123 dexp++; 2124 expCount += dexp; 2125 absValue *= 2.0 ^^ -dexp; 2126 // If the original value was subnormal or if it was a real 2127 // then absValue can still be outside the [0.5, 1.0) range. 2128 if (absValue < 0.5) 2129 { 2130 assert(T.mant_dig > double.mant_dig || isSubnormal(value)); 2131 do 2132 { 2133 absValue += absValue; 2134 expCount--; 2135 } while (absValue < 0.5); 2136 } 2137 else 2138 { 2139 assert(absValue < 1 || T.mant_dig > double.mant_dig); 2140 for (; absValue >= 1; absValue *= T(0.5)) 2141 expCount++; 2142 } 2143 exp = expCount; 2144 return value < 0 ? -absValue : absValue; 2145 } 2146 2147 Unqual!T vf = value; 2148 ushort* vu = cast(ushort*)&vf; 2149 static if (is(immutable T == immutable float)) 2150 int* vi = cast(int*)&vf; 2151 else 2152 long* vl = cast(long*)&vf; 2153 int ex; 2154 alias F = floatTraits!T; 2155 2156 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2157 static if (F.realFormat == RealFormat.ieeeExtended || 2158 F.realFormat == RealFormat.ieeeExtended53) 2159 { 2160 if (ex) 2161 { // If exponent is non-zero 2162 if (ex == F.EXPMASK) // infinity or NaN 2163 { 2164 if (*vl & 0x7FFF_FFFF_FFFF_FFFF) // NaN 2165 { 2166 *vl |= 0xC000_0000_0000_0000; // convert NaNS to NaNQ 2167 exp = int.min; 2168 } 2169 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity 2170 exp = int.min; 2171 else // positive infinity 2172 exp = int.max; 2173 2174 } 2175 else 2176 { 2177 exp = ex - F.EXPBIAS; 2178 vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE; 2179 } 2180 } 2181 else if (!*vl) 2182 { 2183 // vf is +-0.0 2184 exp = 0; 2185 } 2186 else 2187 { 2188 // subnormal 2189 2190 vf *= F.RECIP_EPSILON; 2191 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2192 exp = ex - F.EXPBIAS - T.mant_dig + 1; 2193 vu[F.EXPPOS_SHORT] = ((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FFE; 2194 } 2195 return vf; 2196 } 2197 else static if (F.realFormat == RealFormat.ieeeQuadruple) 2198 { 2199 if (ex) // If exponent is non-zero 2200 { 2201 if (ex == F.EXPMASK) 2202 { 2203 // infinity or NaN 2204 if (vl[MANTISSA_LSB] | 2205 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN 2206 { 2207 // convert NaNS to NaNQ 2208 vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000; 2209 exp = int.min; 2210 } 2211 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity 2212 exp = int.min; 2213 else // positive infinity 2214 exp = int.max; 2215 } 2216 else 2217 { 2218 exp = ex - F.EXPBIAS; 2219 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]); 2220 } 2221 } 2222 else if ((vl[MANTISSA_LSB] | 2223 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) 2224 { 2225 // vf is +-0.0 2226 exp = 0; 2227 } 2228 else 2229 { 2230 // subnormal 2231 vf *= F.RECIP_EPSILON; 2232 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2233 exp = ex - F.EXPBIAS - T.mant_dig + 1; 2234 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]); 2235 } 2236 return vf; 2237 } 2238 else static if (F.realFormat == RealFormat.ieeeDouble) 2239 { 2240 if (ex) // If exponent is non-zero 2241 { 2242 if (ex == F.EXPMASK) // infinity or NaN 2243 { 2244 if (*vl == 0x7FF0_0000_0000_0000) // positive infinity 2245 { 2246 exp = int.max; 2247 } 2248 else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity 2249 exp = int.min; 2250 else 2251 { // NaN 2252 *vl |= 0x0008_0000_0000_0000; // convert NaNS to NaNQ 2253 exp = int.min; 2254 } 2255 } 2256 else 2257 { 2258 exp = (ex - F.EXPBIAS) >> 4; 2259 vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0); 2260 } 2261 } 2262 else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF)) 2263 { 2264 // vf is +-0.0 2265 exp = 0; 2266 } 2267 else 2268 { 2269 // subnormal 2270 vf *= F.RECIP_EPSILON; 2271 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2272 exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1; 2273 vu[F.EXPPOS_SHORT] = 2274 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FE0); 2275 } 2276 return vf; 2277 } 2278 else static if (F.realFormat == RealFormat.ieeeSingle) 2279 { 2280 if (ex) // If exponent is non-zero 2281 { 2282 if (ex == F.EXPMASK) // infinity or NaN 2283 { 2284 if (*vi == 0x7F80_0000) // positive infinity 2285 { 2286 exp = int.max; 2287 } 2288 else if (*vi == 0xFF80_0000) // negative infinity 2289 exp = int.min; 2290 else 2291 { // NaN 2292 *vi |= 0x0040_0000; // convert NaNS to NaNQ 2293 exp = int.min; 2294 } 2295 } 2296 else 2297 { 2298 exp = (ex - F.EXPBIAS) >> 7; 2299 vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00); 2300 } 2301 } 2302 else if (!(*vi & 0x7FFF_FFFF)) 2303 { 2304 // vf is +-0.0 2305 exp = 0; 2306 } 2307 else 2308 { 2309 // subnormal 2310 vf *= F.RECIP_EPSILON; 2311 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2312 exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1; 2313 vu[F.EXPPOS_SHORT] = 2314 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00); 2315 } 2316 return vf; 2317 } 2318 else // static if (F.realFormat == RealFormat.ibmExtended) 2319 { 2320 assert(0, "frexp not implemented"); 2321 } 2322 } 2323 2324 /// 2325 @safe unittest 2326 { 2327 import std.math.operations : isClose; 2328 2329 int exp; 2330 real mantissa = frexp(123.456L, exp); 2331 2332 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L)); 2333 2334 assert(frexp(-real.nan, exp) && exp == int.min); 2335 assert(frexp(real.nan, exp) && exp == int.min); 2336 assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min); 2337 assert(frexp(real.infinity, exp) == real.infinity && exp == int.max); 2338 assert(frexp(-0.0, exp) == -0.0 && exp == 0); 2339 assert(frexp(0.0, exp) == 0.0 && exp == 0); 2340 } 2341 2342 @safe @nogc nothrow unittest 2343 { 2344 import std.math.operations : isClose; 2345 2346 int exp; 2347 real mantissa = frexp(123.456L, exp); 2348 2349 // check if values are equal to 19 decimal digits of precision 2350 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L, 1e-18)); 2351 } 2352 2353 @safe unittest 2354 { 2355 import std.math : floatTraits, RealFormat; 2356 import std.math.traits : isIdentical; 2357 import std.meta : AliasSeq; 2358 import std.typecons : tuple, Tuple; 2359 2360 static foreach (T; AliasSeq!(real, double, float)) 2361 {{ 2362 Tuple!(T, T, int)[] vals = [ // x,frexp,exp 2363 tuple(T(0.0), T( 0.0 ), 0), 2364 tuple(T(-0.0), T( -0.0), 0), 2365 tuple(T(1.0), T( .5 ), 1), 2366 tuple(T(-1.0), T( -.5 ), 1), 2367 tuple(T(2.0), T( .5 ), 2), 2368 tuple(T(float.min_normal/2.0f), T(.5), -126), 2369 tuple(T.infinity, T.infinity, int.max), 2370 tuple(-T.infinity, -T.infinity, int.min), 2371 tuple(T.nan, T.nan, int.min), 2372 tuple(-T.nan, -T.nan, int.min), 2373 2374 // https://issues.dlang.org/show_bug.cgi?id=16026: 2375 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2) 2376 ]; 2377 2378 foreach (elem; vals) 2379 { 2380 T x = elem[0]; 2381 T e = elem[1]; 2382 int exp = elem[2]; 2383 int eptr; 2384 T v = frexp(x, eptr); 2385 assert(isIdentical(e, v)); 2386 assert(exp == eptr); 2387 } 2388 2389 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended) 2390 { 2391 static T[3][] extendedvals = [ // x,frexp,exp 2392 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal 2393 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063], 2394 [T.min_normal, .5, -16381], 2395 [T.min_normal/2.0L, .5, -16382] // subnormal 2396 ]; 2397 foreach (elem; extendedvals) 2398 { 2399 T x = elem[0]; 2400 T e = elem[1]; 2401 int exp = cast(int) elem[2]; 2402 int eptr; 2403 T v = frexp(x, eptr); 2404 assert(isIdentical(e, v)); 2405 assert(exp == eptr); 2406 } 2407 } 2408 }} 2409 2410 // CTFE 2411 alias CtfeFrexpResult= Tuple!(real, int); 2412 static CtfeFrexpResult ctfeFrexp(T)(const T value) 2413 { 2414 int exp; 2415 auto significand = frexp(value, exp); 2416 return CtfeFrexpResult(significand, exp); 2417 } 2418 static foreach (T; AliasSeq!(real, double, float)) 2419 {{ 2420 enum Tuple!(T, T, int)[] vals = [ // x,frexp,exp 2421 tuple(T(0.0), T( 0.0 ), 0), 2422 tuple(T(-0.0), T( -0.0), 0), 2423 tuple(T(1.0), T( .5 ), 1), 2424 tuple(T(-1.0), T( -.5 ), 1), 2425 tuple(T(2.0), T( .5 ), 2), 2426 tuple(T(float.min_normal/2.0f), T(.5), -126), 2427 tuple(T.infinity, T.infinity, int.max), 2428 tuple(-T.infinity, -T.infinity, int.min), 2429 tuple(T.nan, T.nan, int.min), 2430 tuple(-T.nan, -T.nan, int.min), 2431 2432 // https://issues.dlang.org/show_bug.cgi?id=16026: 2433 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2) 2434 ]; 2435 2436 static foreach (elem; vals) 2437 { 2438 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], elem[2])); 2439 } 2440 2441 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended) 2442 { 2443 enum T[3][] extendedvals = [ // x,frexp,exp 2444 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal 2445 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063], 2446 [T.min_normal, .5, -16381], 2447 [T.min_normal/2.0L, .5, -16382] // subnormal 2448 ]; 2449 static foreach (elem; extendedvals) 2450 { 2451 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], cast(int) elem[2])); 2452 } 2453 } 2454 }} 2455 } 2456 2457 @safe unittest 2458 { 2459 import std.meta : AliasSeq; 2460 void foo() { 2461 static foreach (T; AliasSeq!(real, double, float)) 2462 {{ 2463 int exp; 2464 const T a = 1; 2465 immutable T b = 2; 2466 auto c = frexp(a, exp); 2467 auto d = frexp(b, exp); 2468 }} 2469 } 2470 } 2471 2472 /****************************************** 2473 * Extracts the exponent of x as a signed integral value. 2474 * 2475 * If x is not a special value, the result is the same as 2476 * $(D cast(int) logb(x)). 2477 * 2478 * $(TABLE_SV 2479 * $(TR $(TH x) $(TH ilogb(x)) $(TH Range error?)) 2480 * $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes)) 2481 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max) $(TD no)) 2482 * $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD no)) 2483 * ) 2484 */ 2485 int ilogb(T)(const T x) @trusted pure nothrow @nogc 2486 if (isFloatingPoint!T) 2487 { 2488 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB; 2489 2490 import core.bitop : bsr; 2491 alias F = floatTraits!T; 2492 2493 union floatBits 2494 { 2495 T rv; 2496 ushort[T.sizeof/2] vu; 2497 uint[T.sizeof/4] vui; 2498 static if (T.sizeof >= 8) 2499 ulong[T.sizeof/8] vul; 2500 } 2501 floatBits y = void; 2502 y.rv = x; 2503 2504 int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK; 2505 static if (F.realFormat == RealFormat.ieeeExtended || 2506 F.realFormat == RealFormat.ieeeExtended53) 2507 { 2508 if (ex) 2509 { 2510 // If exponent is non-zero 2511 if (ex == F.EXPMASK) // infinity or NaN 2512 { 2513 if (y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) // NaN 2514 return FP_ILOGBNAN; 2515 else // +-infinity 2516 return int.max; 2517 } 2518 else 2519 { 2520 return ex - F.EXPBIAS - 1; 2521 } 2522 } 2523 else if (!y.vul[0]) 2524 { 2525 // vf is +-0.0 2526 return FP_ILOGB0; 2527 } 2528 else 2529 { 2530 // subnormal 2531 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]); 2532 } 2533 } 2534 else static if (F.realFormat == RealFormat.ieeeQuadruple) 2535 { 2536 if (ex) // If exponent is non-zero 2537 { 2538 if (ex == F.EXPMASK) 2539 { 2540 // infinity or NaN 2541 if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN 2542 return FP_ILOGBNAN; 2543 else // +- infinity 2544 return int.max; 2545 } 2546 else 2547 { 2548 return ex - F.EXPBIAS - 1; 2549 } 2550 } 2551 else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) 2552 { 2553 // vf is +-0.0 2554 return FP_ILOGB0; 2555 } 2556 else 2557 { 2558 // subnormal 2559 const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF; 2560 const ulong lsb = y.vul[MANTISSA_LSB]; 2561 if (msb) 2562 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64; 2563 else 2564 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb); 2565 } 2566 } 2567 else static if (F.realFormat == RealFormat.ieeeDouble) 2568 { 2569 if (ex) // If exponent is non-zero 2570 { 2571 if (ex == F.EXPMASK) // infinity or NaN 2572 { 2573 if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000) // +- infinity 2574 return int.max; 2575 else // NaN 2576 return FP_ILOGBNAN; 2577 } 2578 else 2579 { 2580 return ((ex - F.EXPBIAS) >> 4) - 1; 2581 } 2582 } 2583 else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF)) 2584 { 2585 // vf is +-0.0 2586 return FP_ILOGB0; 2587 } 2588 else 2589 { 2590 // subnormal 2591 enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF; 2592 return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64); 2593 } 2594 } 2595 else static if (F.realFormat == RealFormat.ieeeSingle) 2596 { 2597 if (ex) // If exponent is non-zero 2598 { 2599 if (ex == F.EXPMASK) // infinity or NaN 2600 { 2601 if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000) // +- infinity 2602 return int.max; 2603 else // NaN 2604 return FP_ILOGBNAN; 2605 } 2606 else 2607 { 2608 return ((ex - F.EXPBIAS) >> 7) - 1; 2609 } 2610 } 2611 else if (!(y.vui[0] & 0x7FFF_FFFF)) 2612 { 2613 // vf is +-0.0 2614 return FP_ILOGB0; 2615 } 2616 else 2617 { 2618 // subnormal 2619 const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT; 2620 return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa); 2621 } 2622 } 2623 else // static if (F.realFormat == RealFormat.ibmExtended) 2624 { 2625 assert(0, "ilogb not implemented"); 2626 } 2627 } 2628 /// ditto 2629 int ilogb(T)(const T x) @safe pure nothrow @nogc 2630 if (isIntegral!T && isUnsigned!T) 2631 { 2632 import core.bitop : bsr; 2633 if (x == 0) 2634 return FP_ILOGB0; 2635 else 2636 { 2637 static assert(T.sizeof <= ulong.sizeof, "integer size too large for the current ilogb implementation"); 2638 return bsr(x); 2639 } 2640 } 2641 /// ditto 2642 int ilogb(T)(const T x) @safe pure nothrow @nogc 2643 if (isIntegral!T && isSigned!T) 2644 { 2645 import std.traits : Unsigned; 2646 // Note: abs(x) can not be used because the return type is not Unsigned and 2647 // the return value would be wrong for x == int.min 2648 Unsigned!T absx = x >= 0 ? x : -x; 2649 return ilogb(absx); 2650 } 2651 2652 /// 2653 @safe pure unittest 2654 { 2655 assert(ilogb(1) == 0); 2656 assert(ilogb(3) == 1); 2657 assert(ilogb(3.0) == 1); 2658 assert(ilogb(100_000_000) == 26); 2659 2660 assert(ilogb(0) == FP_ILOGB0); 2661 assert(ilogb(0.0) == FP_ILOGB0); 2662 assert(ilogb(double.nan) == FP_ILOGBNAN); 2663 assert(ilogb(double.infinity) == int.max); 2664 } 2665 2666 /** 2667 Special return values of $(LREF ilogb). 2668 */ 2669 alias FP_ILOGB0 = core.stdc.math.FP_ILOGB0; 2670 /// ditto 2671 alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN; 2672 2673 /// 2674 @safe pure unittest 2675 { 2676 assert(ilogb(0) == FP_ILOGB0); 2677 assert(ilogb(0.0) == FP_ILOGB0); 2678 assert(ilogb(double.nan) == FP_ILOGBNAN); 2679 } 2680 2681 @safe nothrow @nogc unittest 2682 { 2683 import std.math : floatTraits, RealFormat; 2684 import std.math.operations : nextUp; 2685 import std.meta : AliasSeq; 2686 import std.typecons : Tuple; 2687 static foreach (F; AliasSeq!(float, double, real)) 2688 {{ 2689 alias T = Tuple!(F, int); 2690 T[13] vals = // x, ilogb(x) 2691 [ 2692 T( F.nan , FP_ILOGBNAN ), 2693 T( -F.nan , FP_ILOGBNAN ), 2694 T( F.infinity, int.max ), 2695 T( -F.infinity, int.max ), 2696 T( 0.0 , FP_ILOGB0 ), 2697 T( -0.0 , FP_ILOGB0 ), 2698 T( 2.0 , 1 ), 2699 T( 2.0001 , 1 ), 2700 T( 1.9999 , 0 ), 2701 T( 0.5 , -1 ), 2702 T( 123.123 , 6 ), 2703 T( -123.123 , 6 ), 2704 T( 0.123 , -4 ), 2705 ]; 2706 2707 foreach (elem; vals) 2708 { 2709 assert(ilogb(elem[0]) == elem[1]); 2710 } 2711 }} 2712 2713 // min_normal and subnormals 2714 assert(ilogb(-float.min_normal) == -126); 2715 assert(ilogb(nextUp(-float.min_normal)) == -127); 2716 assert(ilogb(nextUp(-float(0.0))) == -149); 2717 assert(ilogb(-double.min_normal) == -1022); 2718 assert(ilogb(nextUp(-double.min_normal)) == -1023); 2719 assert(ilogb(nextUp(-double(0.0))) == -1074); 2720 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended) 2721 { 2722 assert(ilogb(-real.min_normal) == -16382); 2723 assert(ilogb(nextUp(-real.min_normal)) == -16383); 2724 assert(ilogb(nextUp(-real(0.0))) == -16445); 2725 } 2726 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) 2727 { 2728 assert(ilogb(-real.min_normal) == -1022); 2729 assert(ilogb(nextUp(-real.min_normal)) == -1023); 2730 assert(ilogb(nextUp(-real(0.0))) == -1074); 2731 } 2732 2733 // test integer types 2734 assert(ilogb(0) == FP_ILOGB0); 2735 assert(ilogb(int.max) == 30); 2736 assert(ilogb(int.min) == 31); 2737 assert(ilogb(uint.max) == 31); 2738 assert(ilogb(long.max) == 62); 2739 assert(ilogb(long.min) == 63); 2740 assert(ilogb(ulong.max) == 63); 2741 } 2742 2743 /******************************************* 2744 * Compute n * 2$(SUPERSCRIPT exp) 2745 * References: frexp 2746 */ 2747 2748 pragma(inline, true) 2749 real ldexp(real n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); } 2750 ///ditto 2751 pragma(inline, true) 2752 double ldexp(double n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); } 2753 ///ditto 2754 pragma(inline, true) 2755 float ldexp(float n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); } 2756 2757 /// 2758 @nogc @safe pure nothrow unittest 2759 { 2760 import std.meta : AliasSeq; 2761 static foreach (T; AliasSeq!(float, double, real)) 2762 {{ 2763 T r; 2764 2765 r = ldexp(3.0L, 3); 2766 assert(r == 24); 2767 2768 r = ldexp(cast(T) 3.0, cast(int) 3); 2769 assert(r == 24); 2770 2771 T n = 3.0; 2772 int exp = 3; 2773 r = ldexp(n, exp); 2774 assert(r == 24); 2775 }} 2776 } 2777 2778 @safe pure nothrow @nogc unittest 2779 { 2780 import std.math : floatTraits, RealFormat; 2781 2782 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended || 2783 floatTraits!(real).realFormat == RealFormat.ieeeExtended53 || 2784 floatTraits!(real).realFormat == RealFormat.ieeeQuadruple) 2785 { 2786 assert(ldexp(1.0L, -16384) == 0x1p-16384L); 2787 assert(ldexp(1.0L, -16382) == 0x1p-16382L); 2788 int x; 2789 real n = frexp(0x1p-16384L, x); 2790 assert(n == 0.5L); 2791 assert(x==-16383); 2792 assert(ldexp(n, x)==0x1p-16384L); 2793 } 2794 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) 2795 { 2796 assert(ldexp(1.0L, -1024) == 0x1p-1024L); 2797 assert(ldexp(1.0L, -1022) == 0x1p-1022L); 2798 int x; 2799 real n = frexp(0x1p-1024L, x); 2800 assert(n == 0.5L); 2801 assert(x==-1023); 2802 assert(ldexp(n, x)==0x1p-1024L); 2803 } 2804 else static assert(false, "Floating point type real not supported"); 2805 } 2806 2807 /* workaround https://issues.dlang.org/show_bug.cgi?id=14718 2808 float parsing depends on platform strtold 2809 @safe pure nothrow @nogc unittest 2810 { 2811 assert(ldexp(1.0, -1024) == 0x1p-1024); 2812 assert(ldexp(1.0, -1022) == 0x1p-1022); 2813 int x; 2814 double n = frexp(0x1p-1024, x); 2815 assert(n == 0.5); 2816 assert(x==-1023); 2817 assert(ldexp(n, x)==0x1p-1024); 2818 } 2819 2820 @safe pure nothrow @nogc unittest 2821 { 2822 assert(ldexp(1.0f, -128) == 0x1p-128f); 2823 assert(ldexp(1.0f, -126) == 0x1p-126f); 2824 int x; 2825 float n = frexp(0x1p-128f, x); 2826 assert(n == 0.5f); 2827 assert(x==-127); 2828 assert(ldexp(n, x)==0x1p-128f); 2829 } 2830 */ 2831 2832 @safe @nogc nothrow unittest 2833 { 2834 import std.math.operations : isClose; 2835 2836 static real[3][] vals = // value,exp,ldexp 2837 [ 2838 [ 0, 0, 0], 2839 [ 1, 0, 1], 2840 [ -1, 0, -1], 2841 [ 1, 1, 2], 2842 [ 123, 10, 125952], 2843 [ real.max, int.max, real.infinity], 2844 [ real.max, -int.max, 0], 2845 [ real.min_normal, -int.max, 0], 2846 ]; 2847 int i; 2848 2849 for (i = 0; i < vals.length; i++) 2850 { 2851 real x = vals[i][0]; 2852 int exp = cast(int) vals[i][1]; 2853 real z = vals[i][2]; 2854 real l = ldexp(x, exp); 2855 2856 assert(isClose(z, l, 1e-6)); 2857 } 2858 2859 real function(real, int) pldexp = &ldexp; 2860 assert(pldexp != null); 2861 } 2862 2863 private 2864 { 2865 import std.math : floatTraits, RealFormat; 2866 2867 version (INLINE_YL2X) {} else 2868 { 2869 static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple) 2870 { 2871 // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x) 2872 static immutable real[13] logCoeffsP = [ 2873 1.313572404063446165910279910527789794488E4L, 2874 7.771154681358524243729929227226708890930E4L, 2875 2.014652742082537582487669938141683759923E5L, 2876 3.007007295140399532324943111654767187848E5L, 2877 2.854829159639697837788887080758954924001E5L, 2878 1.797628303815655343403735250238293741397E5L, 2879 7.594356839258970405033155585486712125861E4L, 2880 2.128857716871515081352991964243375186031E4L, 2881 3.824952356185897735160588078446136783779E3L, 2882 4.114517881637811823002128927449878962058E2L, 2883 2.321125933898420063925789532045674660756E1L, 2884 4.998469661968096229986658302195402690910E-1L, 2885 1.538612243596254322971797716843006400388E-6L 2886 ]; 2887 static immutable real[13] logCoeffsQ = [ 2888 3.940717212190338497730839731583397586124E4L, 2889 2.626900195321832660448791748036714883242E5L, 2890 7.777690340007566932935753241556479363645E5L, 2891 1.347518538384329112529391120390701166528E6L, 2892 1.514882452993549494932585972882995548426E6L, 2893 1.158019977462989115839826904108208787040E6L, 2894 6.132189329546557743179177159925690841200E5L, 2895 2.248234257620569139969141618556349415120E5L, 2896 5.605842085972455027590989944010492125825E4L, 2897 9.147150349299596453976674231612674085381E3L, 2898 9.104928120962988414618126155557301584078E2L, 2899 4.839208193348159620282142911143429644326E1L, 2900 1.0 2901 ]; 2902 2903 // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2) 2904 // where z = 2(x-1)/(x+1) 2905 static immutable real[6] logCoeffsR = [ 2906 1.418134209872192732479751274970992665513E5L, 2907 -8.977257995689735303686582344659576526998E4L, 2908 2.048819892795278657810231591630928516206E4L, 2909 -2.024301798136027039250415126250455056397E3L, 2910 8.057002716646055371965756206836056074715E1L, 2911 -8.828896441624934385266096344596648080902E-1L 2912 ]; 2913 static immutable real[7] logCoeffsS = [ 2914 1.701761051846631278975701529965589676574E6L, 2915 -1.332535117259762928288745111081235577029E6L, 2916 4.001557694070773974936904547424676279307E5L, 2917 -5.748542087379434595104154610899551484314E4L, 2918 3.998526750980007367835804959888064681098E3L, 2919 -1.186359407982897997337150403816839480438E2L, 2920 1.0 2921 ]; 2922 } 2923 else 2924 { 2925 // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x) 2926 static immutable real[7] logCoeffsP = [ 2927 2.0039553499201281259648E1L, 2928 5.7112963590585538103336E1L, 2929 6.0949667980987787057556E1L, 2930 2.9911919328553073277375E1L, 2931 6.5787325942061044846969E0L, 2932 4.9854102823193375972212E-1L, 2933 4.5270000862445199635215E-5L, 2934 ]; 2935 static immutable real[7] logCoeffsQ = [ 2936 6.0118660497603843919306E1L, 2937 2.1642788614495947685003E2L, 2938 3.0909872225312059774938E2L, 2939 2.2176239823732856465394E2L, 2940 8.3047565967967209469434E1L, 2941 1.5062909083469192043167E1L, 2942 1.0000000000000000000000E0L, 2943 ]; 2944 2945 // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2) 2946 // where z = 2(x-1)/(x+1) 2947 static immutable real[4] logCoeffsR = [ 2948 -3.5717684488096787370998E1L, 2949 1.0777257190312272158094E1L, 2950 -7.1990767473014147232598E-1L, 2951 1.9757429581415468984296E-3L, 2952 ]; 2953 static immutable real[4] logCoeffsS = [ 2954 -4.2861221385716144629696E2L, 2955 1.9361891836232102174846E2L, 2956 -2.6201045551331104417768E1L, 2957 1.0000000000000000000000E0L, 2958 ]; 2959 } 2960 } 2961 } 2962 2963 /************************************** 2964 * Calculate the natural logarithm of x. 2965 * 2966 * $(TABLE_SV 2967 * $(TR $(TH x) $(TH log(x)) $(TH divide by 0?) $(TH invalid?)) 2968 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) 2969 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes)) 2970 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) 2971 * ) 2972 */ 2973 real log(real x) @safe pure nothrow @nogc 2974 { 2975 import std.math.constants : LN2, LOG2, SQRT1_2; 2976 import std.math.traits : isInfinity, isNaN, signbit; 2977 import std.math.algebraic : poly; 2978 2979 version (INLINE_YL2X) 2980 return core.math.yl2x(x, LN2); 2981 else 2982 { 2983 // C1 + C2 = LN2. 2984 enum real C1 = 6.93145751953125E-1L; 2985 enum real C2 = 1.428606820309417232121458176568075500134E-6L; 2986 2987 // Special cases. 2988 if (isNaN(x)) 2989 return x; 2990 if (isInfinity(x) && !signbit(x)) 2991 return x; 2992 if (x == 0.0) 2993 return -real.infinity; 2994 if (x < 0.0) 2995 return real.nan; 2996 2997 // Separate mantissa from exponent. 2998 // Note, frexp is used so that denormal numbers will be handled properly. 2999 real y, z; 3000 int exp; 3001 3002 x = frexp(x, exp); 3003 3004 // Logarithm using log(x) = z + z^^3 R(z) / S(z), 3005 // where z = 2(x - 1)/(x + 1) 3006 if ((exp > 2) || (exp < -2)) 3007 { 3008 if (x < SQRT1_2) 3009 { // 2(2x - 1)/(2x + 1) 3010 exp -= 1; 3011 z = x - 0.5; 3012 y = 0.5 * z + 0.5; 3013 } 3014 else 3015 { // 2(x - 1)/(x + 1) 3016 z = x - 0.5; 3017 z -= 0.5; 3018 y = 0.5 * x + 0.5; 3019 } 3020 x = z / y; 3021 z = x * x; 3022 z = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS)); 3023 z += exp * C2; 3024 z += x; 3025 z += exp * C1; 3026 3027 return z; 3028 } 3029 3030 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x) 3031 if (x < SQRT1_2) 3032 { 3033 exp -= 1; 3034 x = 2.0 * x - 1.0; 3035 } 3036 else 3037 { 3038 x = x - 1.0; 3039 } 3040 z = x * x; 3041 y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ)); 3042 y += exp * C2; 3043 z = y - 0.5 * z; 3044 3045 // Note, the sum of above terms does not exceed x/4, 3046 // so it contributes at most about 1/4 lsb to the error. 3047 z += x; 3048 z += exp * C1; 3049 3050 return z; 3051 } 3052 } 3053 3054 /// 3055 @safe pure nothrow @nogc unittest 3056 { 3057 import std.math.operations : feqrel; 3058 import std.math.constants : E; 3059 3060 assert(feqrel(log(E), 1) >= real.mant_dig - 1); 3061 } 3062 3063 /************************************** 3064 * Calculate the base-10 logarithm of x. 3065 * 3066 * $(TABLE_SV 3067 * $(TR $(TH x) $(TH log10(x)) $(TH divide by 0?) $(TH invalid?)) 3068 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) 3069 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes)) 3070 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) 3071 * ) 3072 */ 3073 real log10(real x) @safe pure nothrow @nogc 3074 { 3075 import std.math.constants : LOG2, LN2, SQRT1_2; 3076 import std.math.algebraic : poly; 3077 import std.math.traits : isNaN, isInfinity, signbit; 3078 3079 version (INLINE_YL2X) 3080 return core.math.yl2x(x, LOG2); 3081 else 3082 { 3083 // log10(2) split into two parts. 3084 enum real L102A = 0.3125L; 3085 enum real L102B = -1.14700043360188047862611052755069732318101185E-2L; 3086 3087 // log10(e) split into two parts. 3088 enum real L10EA = 0.5L; 3089 enum real L10EB = -6.570551809674817234887108108339491770560299E-2L; 3090 3091 // Special cases are the same as for log. 3092 if (isNaN(x)) 3093 return x; 3094 if (isInfinity(x) && !signbit(x)) 3095 return x; 3096 if (x == 0.0) 3097 return -real.infinity; 3098 if (x < 0.0) 3099 return real.nan; 3100 3101 // Separate mantissa from exponent. 3102 // Note, frexp is used so that denormal numbers will be handled properly. 3103 real y, z; 3104 int exp; 3105 3106 x = frexp(x, exp); 3107 3108 // Logarithm using log(x) = z + z^^3 R(z) / S(z), 3109 // where z = 2(x - 1)/(x + 1) 3110 if ((exp > 2) || (exp < -2)) 3111 { 3112 if (x < SQRT1_2) 3113 { // 2(2x - 1)/(2x + 1) 3114 exp -= 1; 3115 z = x - 0.5; 3116 y = 0.5 * z + 0.5; 3117 } 3118 else 3119 { // 2(x - 1)/(x + 1) 3120 z = x - 0.5; 3121 z -= 0.5; 3122 y = 0.5 * x + 0.5; 3123 } 3124 x = z / y; 3125 z = x * x; 3126 y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS)); 3127 goto Ldone; 3128 } 3129 3130 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x) 3131 if (x < SQRT1_2) 3132 { 3133 exp -= 1; 3134 x = 2.0 * x - 1.0; 3135 } 3136 else 3137 x = x - 1.0; 3138 3139 z = x * x; 3140 y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ)); 3141 y = y - 0.5 * z; 3142 3143 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2). 3144 // This sequence of operations is critical and it may be horribly 3145 // defeated by some compiler optimizers. 3146 Ldone: 3147 z = y * L10EB; 3148 z += x * L10EB; 3149 z += exp * L102B; 3150 z += y * L10EA; 3151 z += x * L10EA; 3152 z += exp * L102A; 3153 3154 return z; 3155 } 3156 } 3157 3158 /// 3159 @safe pure nothrow @nogc unittest 3160 { 3161 import std.math.algebraic : fabs; 3162 3163 assert(fabs(log10(1000) - 3) < .000001); 3164 } 3165 3166 /** 3167 * Calculates the natural logarithm of 1 + x. 3168 * 3169 * For very small x, log1p(x) will be more accurate than 3170 * log(1 + x). 3171 * 3172 * $(TABLE_SV 3173 * $(TR $(TH x) $(TH log1p(x)) $(TH divide by 0?) $(TH invalid?)) 3174 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) $(TD no)) 3175 * $(TR $(TD -1.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) 3176 * $(TR $(TD $(LT)-1.0) $(TD -$(NAN)) $(TD no) $(TD yes)) 3177 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) 3178 * ) 3179 */ 3180 real log1p(real x) @safe pure nothrow @nogc 3181 { 3182 import std.math.traits : isNaN, isInfinity, signbit; 3183 import std.math.constants : LN2; 3184 3185 version (INLINE_YL2X) 3186 { 3187 // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5, 3188 // ie if -0.29 <= x <= 0.414 3189 return (core.math.fabs(x) <= 0.25) ? core.math.yl2xp1(x, LN2) : core.math.yl2x(x+1, LN2); 3190 } 3191 else 3192 { 3193 // Special cases. 3194 if (isNaN(x) || x == 0.0) 3195 return x; 3196 if (isInfinity(x) && !signbit(x)) 3197 return x; 3198 if (x == -1.0) 3199 return -real.infinity; 3200 if (x < -1.0) 3201 return real.nan; 3202 3203 return log(x + 1.0); 3204 } 3205 } 3206 3207 /// 3208 @safe pure unittest 3209 { 3210 import std.math.traits : isIdentical, isNaN; 3211 import std.math.operations : feqrel; 3212 3213 assert(isIdentical(log1p(0.0), 0.0)); 3214 assert(log1p(1.0).feqrel(0.69314) > 16); 3215 3216 assert(log1p(-1.0) == -real.infinity); 3217 assert(isNaN(log1p(-2.0))); 3218 assert(log1p(real.nan) is real.nan); 3219 assert(log1p(-real.nan) is -real.nan); 3220 assert(log1p(real.infinity) == real.infinity); 3221 } 3222 3223 /*************************************** 3224 * Calculates the base-2 logarithm of x: 3225 * $(SUB log, 2)x 3226 * 3227 * $(TABLE_SV 3228 * $(TR $(TH x) $(TH log2(x)) $(TH divide by 0?) $(TH invalid?)) 3229 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no) ) 3230 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes) ) 3231 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no) ) 3232 * ) 3233 */ 3234 real log2(real x) @safe pure nothrow @nogc 3235 { 3236 import std.math.traits : isNaN, isInfinity, signbit; 3237 import std.math.constants : SQRT1_2, LOG2E; 3238 import std.math.algebraic : poly; 3239 3240 version (INLINE_YL2X) 3241 return core.math.yl2x(x, 1.0L); 3242 else 3243 { 3244 // Special cases are the same as for log. 3245 if (isNaN(x)) 3246 return x; 3247 if (isInfinity(x) && !signbit(x)) 3248 return x; 3249 if (x == 0.0) 3250 return -real.infinity; 3251 if (x < 0.0) 3252 return real.nan; 3253 3254 // Separate mantissa from exponent. 3255 // Note, frexp is used so that denormal numbers will be handled properly. 3256 real y, z; 3257 int exp; 3258 3259 x = frexp(x, exp); 3260 3261 // Logarithm using log(x) = z + z^^3 R(z) / S(z), 3262 // where z = 2(x - 1)/(x + 1) 3263 if ((exp > 2) || (exp < -2)) 3264 { 3265 if (x < SQRT1_2) 3266 { // 2(2x - 1)/(2x + 1) 3267 exp -= 1; 3268 z = x - 0.5; 3269 y = 0.5 * z + 0.5; 3270 } 3271 else 3272 { // 2(x - 1)/(x + 1) 3273 z = x - 0.5; 3274 z -= 0.5; 3275 y = 0.5 * x + 0.5; 3276 } 3277 x = z / y; 3278 z = x * x; 3279 y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS)); 3280 goto Ldone; 3281 } 3282 3283 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x) 3284 if (x < SQRT1_2) 3285 { 3286 exp -= 1; 3287 x = 2.0 * x - 1.0; 3288 } 3289 else 3290 x = x - 1.0; 3291 3292 z = x * x; 3293 y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ)); 3294 y = y - 0.5 * z; 3295 3296 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2). 3297 // This sequence of operations is critical and it may be horribly 3298 // defeated by some compiler optimizers. 3299 Ldone: 3300 z = y * (LOG2E - 1.0); 3301 z += x * (LOG2E - 1.0); 3302 z += y; 3303 z += x; 3304 z += exp; 3305 3306 return z; 3307 } 3308 } 3309 3310 /// 3311 @safe unittest 3312 { 3313 import std.math.operations : isClose; 3314 3315 assert(isClose(log2(1024.0L), 10)); 3316 } 3317 3318 @safe @nogc nothrow unittest 3319 { 3320 import std.math.operations : isClose; 3321 3322 // check if values are equal to 19 decimal digits of precision 3323 assert(isClose(log2(1024.0L), 10, 1e-18)); 3324 } 3325 3326 /***************************************** 3327 * Extracts the exponent of x as a signed integral value. 3328 * 3329 * If x is subnormal, it is treated as if it were normalized. 3330 * For a positive, finite x: 3331 * 3332 * 1 $(LT)= $(I x) * FLT_RADIX$(SUPERSCRIPT -logb(x)) $(LT) FLT_RADIX 3333 * 3334 * $(TABLE_SV 3335 * $(TR $(TH x) $(TH logb(x)) $(TH divide by 0?) ) 3336 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no)) 3337 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) ) 3338 * ) 3339 */ 3340 real logb(real x) @trusted nothrow @nogc 3341 { 3342 version (InlineAsm_X87_MSVC) 3343 { 3344 version (X86_64) 3345 { 3346 asm pure nothrow @nogc 3347 { 3348 naked ; 3349 fld real ptr [RCX] ; 3350 fxtract ; 3351 fstp ST(0) ; 3352 ret ; 3353 } 3354 } 3355 else 3356 { 3357 asm pure nothrow @nogc 3358 { 3359 fld x ; 3360 fxtract ; 3361 fstp ST(0) ; 3362 } 3363 } 3364 } 3365 else 3366 return core.stdc.math.logbl(x); 3367 } 3368 3369 /// 3370 @safe @nogc nothrow unittest 3371 { 3372 assert(logb(1.0) == 0); 3373 assert(logb(100.0) == 6); 3374 3375 assert(logb(0.0) == -real.infinity); 3376 assert(logb(real.infinity) == real.infinity); 3377 assert(logb(-real.infinity) == real.infinity); 3378 } 3379 3380 /************************************* 3381 * Efficiently calculates x * 2$(SUPERSCRIPT n). 3382 * 3383 * scalbn handles underflow and overflow in 3384 * the same fashion as the basic arithmetic operators. 3385 * 3386 * $(TABLE_SV 3387 * $(TR $(TH x) $(TH scalb(x))) 3388 * $(TR $(TD $(PLUSMNINF)) $(TD $(PLUSMNINF)) ) 3389 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) ) 3390 * ) 3391 */ 3392 pragma(inline, true) 3393 real scalbn(real x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); } 3394 3395 /// ditto 3396 pragma(inline, true) 3397 double scalbn(double x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); } 3398 3399 /// ditto 3400 pragma(inline, true) 3401 float scalbn(float x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); } 3402 3403 /// 3404 @safe pure nothrow @nogc unittest 3405 { 3406 assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L); 3407 assert(scalbn(-real.infinity, 5) == -real.infinity); 3408 assert(scalbn(2.0,10) == 2048.0); 3409 assert(scalbn(2048.0f,-10) == 2.0f); 3410 } 3411 3412 pragma(inline, true) 3413 private F _scalbn(F)(F x, int n) 3414 { 3415 import std.math.traits : isInfinity; 3416 3417 if (__ctfe) 3418 { 3419 // Handle special cases. 3420 if (x == F(0.0) || isInfinity(x)) 3421 return x; 3422 } 3423 return core.math.ldexp(x, n); 3424 } 3425 3426 @safe pure nothrow @nogc unittest 3427 { 3428 // CTFE-able test 3429 static assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L); 3430 static assert(scalbn(-real.infinity, 5) == -real.infinity); 3431 // Test with large exponent delta n where the result is in bounds but 2.0L ^^ n is not. 3432 enum initialExponent = real.min_exp + 2, resultExponent = real.max_exp - 2; 3433 enum int n = resultExponent - initialExponent; 3434 enum real x = 0x1.2345678abcdefp0L * (2.0L ^^ initialExponent); 3435 enum staticResult = scalbn(x, n); 3436 static assert(staticResult == 0x1.2345678abcdefp0L * (2.0L ^^ resultExponent)); 3437 assert(scalbn(x, n) == staticResult); 3438 } 3439 3440
Indexes created Tue Feb 24 08:35:24 UTC 2026