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@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 (OSX) { } // macOS 13 (M1) has issues emulating instruction 43 else version = INLINE_YL2X; // x87 has opcodes for these 44 } 45 46 version (D_InlineAsm_X86) version = InlineAsm_X86_Any; 47 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any; 48 49 version (InlineAsm_X86_Any) version = InlineAsm_X87; 50 version (InlineAsm_X87) 51 { 52 static assert(real.mant_dig == 64); 53 version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC; 54 } 55 56 version (D_HardFloat) 57 { 58 // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport 59 version (IeeeFlagsSupport) version = FloatingPointControlSupport; 60 } 61 62 /** 63 * Compute the value of x $(SUPERSCRIPT n), where n is an integer 64 */ 65 Unqual!F pow(F, G)(F x, G n) @nogc @trusted pure nothrow 66 if (isFloatingPoint!(F) && isIntegral!(G)) 67 { 68 import std.traits : Unsigned; 69 70 real p = 1.0, v = void; 71 Unsigned!(Unqual!G) m = n; 72 73 if (n < 0) 74 { 75 if (n == -1) return 1 / x; 76 77 m = cast(typeof(m))(0 - n); 78 v = p / x; 79 } 80 else 81 { 82 switch (n) 83 { 84 case 0: 85 return 1.0; 86 case 1: 87 return x; 88 case 2: 89 return x * x; 90 default: 91 } 92 93 v = x; 94 } 95 96 while (1) 97 { 98 if (m & 1) 99 p *= v; 100 m >>= 1; 101 if (!m) 102 break; 103 v *= v; 104 } 105 return p; 106 } 107 108 /// 109 @safe pure nothrow @nogc unittest 110 { 111 import std.math.operations : feqrel; 112 113 assert(pow(2.0, 5) == 32.0); 114 assert(pow(1.5, 9).feqrel(38.4433) > 16); 115 assert(pow(real.nan, 2) is real.nan); 116 assert(pow(real.infinity, 2) == real.infinity); 117 } 118 119 @safe pure nothrow @nogc unittest 120 { 121 import std.math.operations : isClose, feqrel; 122 123 // Make sure it instantiates and works properly on immutable values and 124 // with various integer and float types. 125 immutable real x = 46; 126 immutable float xf = x; 127 immutable double xd = x; 128 immutable uint one = 1; 129 immutable ushort two = 2; 130 immutable ubyte three = 3; 131 immutable ulong eight = 8; 132 133 immutable int neg1 = -1; 134 immutable short neg2 = -2; 135 immutable byte neg3 = -3; 136 immutable long neg8 = -8; 137 138 139 assert(pow(x,0) == 1.0); 140 assert(pow(xd,one) == x); 141 assert(pow(xf,two) == x * x); 142 assert(pow(x,three) == x * x * x); 143 assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x)); 144 145 assert(pow(x, neg1) == 1 / x); 146 147 assert(isClose(pow(xd, neg2), cast(double) (1 / (x * x)), 1e-15)); 148 assert(isClose(pow(xf, neg8), cast(float) (1 / ((x * x) * (x * x) * (x * x) * (x * x))), 1e-15)); 149 150 assert(feqrel(pow(x, neg3), 1 / (x * x * x)) >= real.mant_dig - 1); 151 } 152 153 @safe @nogc nothrow unittest 154 { 155 import std.math.operations : isClose; 156 157 assert(isClose(pow(2.0L, 10L), 1024, 1e-18)); 158 } 159 160 // https://issues.dlang.org/show_bug.cgi?id=21601 161 @safe @nogc nothrow pure unittest 162 { 163 // When reals are large enough the results of pow(b, e) can be 164 // calculated correctly, if b is of type float or double and e is 165 // not too large. 166 static if (real.mant_dig >= 64) 167 { 168 // expected result: 3.790e-42 169 assert(pow(-513645318757045764096.0f, -2) > 0.0); 170 171 // expected result: 3.763915357831797e-309 172 assert(pow(-1.6299717435255677e+154, -2) > 0.0); 173 } 174 } 175 176 @safe @nogc nothrow unittest 177 { 178 import std.math.operations : isClose; 179 import std.math.traits : isInfinity; 180 181 static float f1 = 19100.0f; 182 static float f2 = 0.000012f; 183 184 assert(isClose(pow(f1,9), 3.3829868e+38f)); 185 assert(isInfinity(pow(f1,10))); 186 assert(pow(f2,9) > 0.0f); 187 assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal)); 188 189 static double d1 = 21800.0; 190 static double d2 = 0.000012; 191 192 assert(isClose(pow(d1,71), 1.0725339442974e+308)); 193 assert(isInfinity(pow(d1,72))); 194 assert(pow(d2,65) > 0.0f); 195 assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal)); 196 197 static if (real.mant_dig == 64) // x87 198 { 199 static real r1 = 21950.0L; 200 static real r2 = 0.000011L; 201 202 assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L)); 203 assert(isInfinity(pow(r1,1137))); 204 assert(pow(r2,998) > 0.0L); 205 assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal)); 206 } 207 } 208 209 @safe @nogc nothrow pure unittest 210 { 211 import std.math.operations : isClose; 212 213 enum f1 = 19100.0f; 214 enum f2 = 0.000012f; 215 216 static assert(isClose(pow(f1,9), 3.3829868e+38f)); 217 static assert(pow(f1,10) > float.max); 218 static assert(pow(f2,9) > 0.0f); 219 static assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal)); 220 221 enum d1 = 21800.0; 222 enum d2 = 0.000012; 223 224 static assert(isClose(pow(d1,71), 1.0725339442974e+308)); 225 static assert(pow(d1,72) > double.max); 226 static assert(pow(d2,65) > 0.0f); 227 static assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal)); 228 229 static if (real.mant_dig == 64) // x87 230 { 231 enum r1 = 21950.0L; 232 enum r2 = 0.000011L; 233 234 static assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L)); 235 static assert(pow(r1,1137) > real.max); 236 static assert(pow(r2,998) > 0.0L); 237 static assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal)); 238 } 239 } 240 241 /** 242 * Compute the power of two integral numbers. 243 * 244 * Params: 245 * x = base 246 * n = exponent 247 * 248 * Returns: 249 * x raised to the power of n. If n is negative the result is 1 / pow(x, -n), 250 * which is calculated as integer division with remainder. This may result in 251 * a division by zero error. 252 * 253 * If both x and n are 0, the result is 1. 254 * 255 * Throws: 256 * If x is 0 and n is negative, the result is the same as the result of a 257 * division by zero. 258 */ 259 typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @nogc @trusted pure nothrow 260 if (isIntegral!(F) && isIntegral!(G)) 261 { 262 import std.traits : isSigned; 263 264 typeof(return) p, v = void; 265 Unqual!G m = n; 266 267 static if (isSigned!(F)) 268 { 269 if (x == -1) return cast(typeof(return)) (m & 1 ? -1 : 1); 270 } 271 static if (isSigned!(G)) 272 { 273 if (x == 0 && m <= -1) return x / 0; 274 } 275 if (x == 1) return 1; 276 static if (isSigned!(G)) 277 { 278 if (m < 0) return 0; 279 } 280 281 switch (m) 282 { 283 case 0: 284 p = 1; 285 break; 286 287 case 1: 288 p = x; 289 break; 290 291 case 2: 292 p = x * x; 293 break; 294 295 default: 296 v = x; 297 p = 1; 298 while (1) 299 { 300 if (m & 1) 301 p *= v; 302 m >>= 1; 303 if (!m) 304 break; 305 v *= v; 306 } 307 break; 308 } 309 return p; 310 } 311 312 /// 313 @safe pure nothrow @nogc unittest 314 { 315 assert(pow(2, 3) == 8); 316 assert(pow(3, 2) == 9); 317 318 assert(pow(2, 10) == 1_024); 319 assert(pow(2, 20) == 1_048_576); 320 assert(pow(2, 30) == 1_073_741_824); 321 322 assert(pow(0, 0) == 1); 323 324 assert(pow(1, -5) == 1); 325 assert(pow(1, -6) == 1); 326 assert(pow(-1, -5) == -1); 327 assert(pow(-1, -6) == 1); 328 329 assert(pow(-2, 5) == -32); 330 assert(pow(-2, -5) == 0); 331 assert(pow(cast(double) -2, -5) == -0.03125); 332 } 333 334 @safe pure nothrow @nogc unittest 335 { 336 immutable int one = 1; 337 immutable byte two = 2; 338 immutable ubyte three = 3; 339 immutable short four = 4; 340 immutable long ten = 10; 341 342 assert(pow(two, three) == 8); 343 assert(pow(two, ten) == 1024); 344 assert(pow(one, ten) == 1); 345 assert(pow(ten, four) == 10_000); 346 assert(pow(four, 10) == 1_048_576); 347 assert(pow(three, four) == 81); 348 } 349 350 // https://issues.dlang.org/show_bug.cgi?id=7006 351 @safe pure nothrow @nogc unittest 352 { 353 assert(pow(5, -1) == 0); 354 assert(pow(-5, -1) == 0); 355 assert(pow(5, -2) == 0); 356 assert(pow(-5, -2) == 0); 357 assert(pow(-1, int.min) == 1); 358 assert(pow(-2, int.min) == 0); 359 360 assert(pow(4294967290UL,2) == 18446744022169944100UL); 361 assert(pow(0,uint.max) == 0); 362 } 363 364 /**Computes integer to floating point powers.*/ 365 real pow(I, F)(I x, F y) @nogc @trusted pure nothrow 366 if (isIntegral!I && isFloatingPoint!F) 367 { 368 return pow(cast(real) x, cast(Unqual!F) y); 369 } 370 371 /// 372 @safe pure nothrow @nogc unittest 373 { 374 assert(pow(2, 5.0) == 32.0); 375 assert(pow(7, 3.0) == 343.0); 376 assert(pow(2, real.nan) is real.nan); 377 assert(pow(2, real.infinity) == real.infinity); 378 } 379 380 /** 381 * Calculates x$(SUPERSCRIPT y). 382 * 383 * $(TABLE_SV 384 * $(TR $(TH x) $(TH y) $(TH pow(x, y)) 385 * $(TH div 0) $(TH invalid?)) 386 * $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD 1.0) 387 * $(TD no) $(TD no) ) 388 * $(TR $(TD |x| $(GT) 1) $(TD +$(INFIN)) $(TD +$(INFIN)) 389 * $(TD no) $(TD no) ) 390 * $(TR $(TD |x| $(LT) 1) $(TD +$(INFIN)) $(TD +0.0) 391 * $(TD no) $(TD no) ) 392 * $(TR $(TD |x| $(GT) 1) $(TD -$(INFIN)) $(TD +0.0) 393 * $(TD no) $(TD no) ) 394 * $(TR $(TD |x| $(LT) 1) $(TD -$(INFIN)) $(TD +$(INFIN)) 395 * $(TD no) $(TD no) ) 396 * $(TR $(TD +$(INFIN)) $(TD $(GT) 0.0) $(TD +$(INFIN)) 397 * $(TD no) $(TD no) ) 398 * $(TR $(TD +$(INFIN)) $(TD $(LT) 0.0) $(TD +0.0) 399 * $(TD no) $(TD no) ) 400 * $(TR $(TD -$(INFIN)) $(TD odd integer $(GT) 0.0) $(TD -$(INFIN)) 401 * $(TD no) $(TD no) ) 402 * $(TR $(TD -$(INFIN)) $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN)) 403 * $(TD no) $(TD no)) 404 * $(TR $(TD -$(INFIN)) $(TD odd integer $(LT) 0.0) $(TD -0.0) 405 * $(TD no) $(TD no) ) 406 * $(TR $(TD -$(INFIN)) $(TD $(LT) 0.0, not odd integer) $(TD +0.0) 407 * $(TD no) $(TD no) ) 408 * $(TR $(TD $(PLUSMN)1.0) $(TD $(PLUSMN)$(INFIN)) $(TD -$(NAN)) 409 * $(TD no) $(TD yes) ) 410 * $(TR $(TD $(LT) 0.0) $(TD finite, nonintegral) $(TD $(NAN)) 411 * $(TD no) $(TD yes)) 412 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(LT) 0.0) $(TD $(PLUSMNINF)) 413 * $(TD yes) $(TD no) ) 414 * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN)) 415 * $(TD yes) $(TD no)) 416 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(GT) 0.0) $(TD $(PLUSMN)0.0) 417 * $(TD no) $(TD no) ) 418 * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT) 0.0, not odd integer) $(TD +0.0) 419 * $(TD no) $(TD no) ) 420 * ) 421 */ 422 Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow 423 if (isFloatingPoint!(F) && isFloatingPoint!(G)) 424 { 425 import core.math : fabs, sqrt; 426 import std.math.traits : isInfinity, isNaN, signbit; 427 428 alias Float = typeof(return); 429 430 static real impl(real x, real y) @nogc pure nothrow 431 { 432 // Special cases. 433 if (isNaN(y)) 434 return y; 435 if (isNaN(x) && y != 0.0) 436 return x; 437 438 // Even if x is NaN. 439 if (y == 0.0) 440 return 1.0; 441 if (y == 1.0) 442 return x; 443 444 if (isInfinity(y)) 445 { 446 if (isInfinity(x)) 447 { 448 if (!signbit(y) && !signbit(x)) 449 return F.infinity; 450 else 451 return F.nan; 452 } 453 else if (fabs(x) > 1) 454 { 455 if (signbit(y)) 456 return +0.0; 457 else 458 return F.infinity; 459 } 460 else if (fabs(x) == 1) 461 { 462 return F.nan; 463 } 464 else // < 1 465 { 466 if (signbit(y)) 467 return F.infinity; 468 else 469 return +0.0; 470 } 471 } 472 if (isInfinity(x)) 473 { 474 if (signbit(x)) 475 { 476 long i = cast(long) y; 477 if (y > 0.0) 478 { 479 if (i == y && i & 1) 480 return -F.infinity; 481 else if (i == y) 482 return F.infinity; 483 else 484 return -F.nan; 485 } 486 else if (y < 0.0) 487 { 488 if (i == y && i & 1) 489 return -0.0; 490 else if (i == y) 491 return +0.0; 492 else 493 return F.nan; 494 } 495 } 496 else 497 { 498 if (y > 0.0) 499 return F.infinity; 500 else if (y < 0.0) 501 return +0.0; 502 } 503 } 504 505 if (x == 0.0) 506 { 507 if (signbit(x)) 508 { 509 long i = cast(long) y; 510 if (y > 0.0) 511 { 512 if (i == y && i & 1) 513 return -0.0; 514 else 515 return +0.0; 516 } 517 else if (y < 0.0) 518 { 519 if (i == y && i & 1) 520 return -F.infinity; 521 else 522 return F.infinity; 523 } 524 } 525 else 526 { 527 if (y > 0.0) 528 return +0.0; 529 else if (y < 0.0) 530 return F.infinity; 531 } 532 } 533 if (x == 1.0) 534 return 1.0; 535 536 if (y >= F.max) 537 { 538 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0)) 539 return 0.0; 540 if (x > 1.0 || x < -1.0) 541 return F.infinity; 542 } 543 if (y <= -F.max) 544 { 545 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0)) 546 return F.infinity; 547 if (x > 1.0 || x < -1.0) 548 return 0.0; 549 } 550 551 if (x >= F.max) 552 { 553 if (y > 0.0) 554 return F.infinity; 555 else 556 return 0.0; 557 } 558 if (x <= -F.max) 559 { 560 long i = cast(long) y; 561 if (y > 0.0) 562 { 563 if (i == y && i & 1) 564 return -F.infinity; 565 else 566 return F.infinity; 567 } 568 else if (y < 0.0) 569 { 570 if (i == y && i & 1) 571 return -0.0; 572 else 573 return +0.0; 574 } 575 } 576 577 // Integer power of x. 578 long iy = cast(long) y; 579 if (iy == y && fabs(y) < 32_768.0) 580 return pow(x, iy); 581 582 real sign = 1.0; 583 if (x < 0) 584 { 585 // Result is real only if y is an integer 586 // Check for a non-zero fractional part 587 enum maxOdd = pow(2.0L, real.mant_dig) - 1.0L; 588 static if (maxOdd > ulong.max) 589 { 590 // Generic method, for any FP type 591 import std.math.rounding : floor; 592 if (floor(y) != y) 593 return sqrt(x); // Complex result -- create a NaN 594 595 const hy = 0.5 * y; 596 if (floor(hy) != hy) 597 sign = -1.0; 598 } 599 else 600 { 601 // Much faster, if ulong has enough precision 602 const absY = fabs(y); 603 if (absY <= maxOdd) 604 { 605 const uy = cast(ulong) absY; 606 if (uy != absY) 607 return sqrt(x); // Complex result -- create a NaN 608 609 if (uy & 1) 610 sign = -1.0; 611 } 612 } 613 x = -x; 614 } 615 version (INLINE_YL2X) 616 { 617 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) ) 618 // TODO: This is not accurate in practice. A fast and accurate 619 // (though complicated) method is described in: 620 // "An efficient rounding boundary test for pow(x, y) 621 // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007). 622 return sign * exp2( core.math.yl2x(x, y) ); 623 } 624 else 625 { 626 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) ) 627 // TODO: This is not accurate in practice. A fast and accurate 628 // (though complicated) method is described in: 629 // "An efficient rounding boundary test for pow(x, y) 630 // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007). 631 Float w = exp2(y * log2(x)); 632 return sign * w; 633 } 634 } 635 return impl(x, y); 636 } 637 638 /// 639 @safe pure nothrow @nogc unittest 640 { 641 import std.math.operations : isClose; 642 643 assert(isClose(pow(2.0, 3.0), 8.0)); 644 assert(isClose(pow(1.5, 10.0), 57.6650390625)); 645 646 // square root of 9 647 assert(isClose(pow(9.0, 0.5), 3.0)); 648 // 10th root of 1024 649 assert(isClose(pow(1024.0, 0.1), 2.0)); 650 651 assert(isClose(pow(-4.0, 3.0), -64.0)); 652 653 // reciprocal of 4 ^^ 2 654 assert(isClose(pow(4.0, -2.0), 0.0625)); 655 // reciprocal of (-2) ^^ 3 656 assert(isClose(pow(-2.0, -3.0), -0.125)); 657 658 assert(isClose(pow(-2.5, 3.0), -15.625)); 659 // reciprocal of 2.5 ^^ 3 660 assert(isClose(pow(2.5, -3.0), 0.064)); 661 // reciprocal of (-2.5) ^^ 3 662 assert(isClose(pow(-2.5, -3.0), -0.064)); 663 664 // reciprocal of square root of 4 665 assert(isClose(pow(4.0, -0.5), 0.5)); 666 667 // per definition 668 assert(isClose(pow(0.0, 0.0), 1.0)); 669 } 670 671 /// 672 @safe pure nothrow @nogc unittest 673 { 674 import std.math.operations : isClose; 675 676 // the result is a complex number 677 // which cannot be represented as floating point number 678 import std.math.traits : isNaN; 679 assert(isNaN(pow(-2.5, -1.5))); 680 681 // use the ^^-operator of std.complex instead 682 import std.complex : complex; 683 auto c1 = complex(-2.5, 0.0); 684 auto c2 = complex(-1.5, 0.0); 685 auto result = c1 ^^ c2; 686 // exact result apparently depends on `real` precision => increased tolerance 687 assert(isClose(result.re, -4.64705438e-17, 2e-4)); 688 assert(isClose(result.im, 2.52982e-1, 2e-4)); 689 } 690 691 @safe pure nothrow @nogc unittest 692 { 693 import std.math.traits : isNaN; 694 695 assert(pow(1.5, real.infinity) == real.infinity); 696 assert(pow(0.5, real.infinity) == 0.0); 697 assert(pow(1.5, -real.infinity) == 0.0); 698 assert(pow(0.5, -real.infinity) == real.infinity); 699 assert(pow(real.infinity, 1.0) == real.infinity); 700 assert(pow(real.infinity, -1.0) == 0.0); 701 assert(pow(real.infinity, real.infinity) == real.infinity); 702 assert(pow(-real.infinity, 1.0) == -real.infinity); 703 assert(pow(-real.infinity, 2.0) == real.infinity); 704 assert(pow(-real.infinity, -1.0) == -0.0); 705 assert(pow(-real.infinity, -2.0) == 0.0); 706 assert(isNaN(pow(1.0, real.infinity))); 707 assert(pow(0.0, -1.0) == real.infinity); 708 assert(pow(real.nan, 0.0) == 1.0); 709 assert(isNaN(pow(real.nan, 3.0))); 710 assert(isNaN(pow(3.0, real.nan))); 711 } 712 713 @safe @nogc nothrow unittest 714 { 715 import std.math.operations : isClose; 716 717 assert(isClose(pow(2.0L, 10.0L), 1024, 1e-18)); 718 } 719 720 @safe pure nothrow @nogc unittest 721 { 722 import std.math.operations : isClose; 723 import std.math.traits : isIdentical, isNaN; 724 import std.math.constants : PI; 725 726 // Test all the special values. These unittests can be run on Windows 727 // by temporarily changing the version (linux) to version (all). 728 immutable float zero = 0; 729 immutable real one = 1; 730 immutable double two = 2; 731 immutable float three = 3; 732 immutable float fnan = float.nan; 733 immutable double dnan = double.nan; 734 immutable real rnan = real.nan; 735 immutable dinf = double.infinity; 736 immutable rninf = -real.infinity; 737 738 assert(pow(fnan, zero) == 1); 739 assert(pow(dnan, zero) == 1); 740 assert(pow(rnan, zero) == 1); 741 742 assert(pow(two, dinf) == double.infinity); 743 assert(isIdentical(pow(0.2f, dinf), +0.0)); 744 assert(pow(0.99999999L, rninf) == real.infinity); 745 assert(isIdentical(pow(1.000000001, rninf), +0.0)); 746 assert(pow(dinf, 0.001) == dinf); 747 assert(isIdentical(pow(dinf, -0.001), +0.0)); 748 assert(pow(rninf, 3.0L) == rninf); 749 assert(pow(rninf, 2.0L) == real.infinity); 750 assert(isIdentical(pow(rninf, -3.0), -0.0)); 751 assert(isIdentical(pow(rninf, -2.0), +0.0)); 752 753 // @@@BUG@@@ somewhere 754 version (OSX) {} else assert(isNaN(pow(one, dinf))); 755 version (OSX) {} else assert(isNaN(pow(-one, dinf))); 756 assert(isNaN(pow(-0.2, PI))); 757 // boundary cases. Note that epsilon == 2^^-n for some n, 758 // so 1/epsilon == 2^^n is always even. 759 assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L); 760 assert(pow(-1.0L, 1/real.epsilon) == 1.0L); 761 assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L))); 762 assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L))); 763 764 assert(pow(0.0, -3.0) == double.infinity); 765 assert(pow(-0.0, -3.0) == -double.infinity); 766 assert(pow(0.0, -PI) == double.infinity); 767 assert(pow(-0.0, -PI) == double.infinity); 768 assert(isIdentical(pow(0.0, 5.0), 0.0)); 769 assert(isIdentical(pow(-0.0, 5.0), -0.0)); 770 assert(isIdentical(pow(0.0, 6.0), 0.0)); 771 assert(isIdentical(pow(-0.0, 6.0), 0.0)); 772 773 // https://issues.dlang.org/show_bug.cgi?id=14786 fixed 774 immutable real maxOdd = pow(2.0L, real.mant_dig) - 1.0L; 775 assert(pow(-1.0L, maxOdd) == -1.0L); 776 assert(pow(-1.0L, -maxOdd) == -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 assert(pow(-1.0L, -maxOdd - 1.0L) == 1.0L); 781 782 // Now, actual numbers. 783 assert(isClose(pow(two, three), 8.0)); 784 assert(isClose(pow(two, -2.5), 0.1767766953)); 785 786 // Test integer to float power. 787 immutable uint twoI = 2; 788 assert(isClose(pow(twoI, three), 8.0)); 789 } 790 791 // https://issues.dlang.org/show_bug.cgi?id=20508 792 @safe pure nothrow @nogc unittest 793 { 794 import std.math.traits : isNaN; 795 796 assert(isNaN(pow(-double.infinity, 0.5))); 797 798 assert(isNaN(pow(-real.infinity, real.infinity))); 799 assert(isNaN(pow(-real.infinity, -real.infinity))); 800 assert(isNaN(pow(-real.infinity, 1.234))); 801 assert(isNaN(pow(-real.infinity, -0.751))); 802 assert(pow(-real.infinity, 0.0) == 1.0); 803 } 804 805 /** Computes the value of a positive integer `x`, raised to the power `n`, modulo `m`. 806 * 807 * Params: 808 * x = base 809 * n = exponent 810 * m = modulus 811 * 812 * Returns: 813 * `x` to the power `n`, modulo `m`. 814 * The return type is the largest of `x`'s and `m`'s type. 815 * 816 * The function requires that all values have unsigned types. 817 */ 818 Unqual!(Largest!(F, H)) powmod(F, G, H)(F x, G n, H m) 819 if (isUnsigned!F && isUnsigned!G && isUnsigned!H) 820 { 821 import std.meta : AliasSeq; 822 823 alias T = Unqual!(Largest!(F, H)); 824 static if (T.sizeof <= 4) 825 { 826 alias DoubleT = AliasSeq!(void, ushort, uint, void, ulong)[T.sizeof]; 827 } 828 829 static T mulmod(T a, T b, T c) 830 { 831 static if (T.sizeof == 8) 832 { 833 static T addmod(T a, T b, T c) 834 { 835 b = c - b; 836 if (a >= b) 837 return a - b; 838 else 839 return c - b + a; 840 } 841 842 T result = 0, tmp; 843 844 b %= c; 845 while (a > 0) 846 { 847 if (a & 1) 848 result = addmod(result, b, c); 849 850 a >>= 1; 851 b = addmod(b, b, c); 852 } 853 854 return result; 855 } 856 else 857 { 858 DoubleT result = cast(DoubleT) (cast(DoubleT) a * cast(DoubleT) b); 859 return result % c; 860 } 861 } 862 863 T base = x, result = 1, modulus = m; 864 Unqual!G exponent = n; 865 866 while (exponent > 0) 867 { 868 if (exponent & 1) 869 result = mulmod(result, base, modulus); 870 871 base = mulmod(base, base, modulus); 872 exponent >>= 1; 873 } 874 875 return result; 876 } 877 878 /// 879 @safe pure nothrow @nogc unittest 880 { 881 assert(powmod(1U, 10U, 3U) == 1); 882 assert(powmod(3U, 2U, 6U) == 3); 883 assert(powmod(5U, 5U, 15U) == 5); 884 assert(powmod(2U, 3U, 5U) == 3); 885 assert(powmod(2U, 4U, 5U) == 1); 886 assert(powmod(2U, 5U, 5U) == 2); 887 } 888 889 @safe pure nothrow @nogc unittest 890 { 891 ulong a = 18446744073709551615u, b = 20u, c = 18446744073709551610u; 892 assert(powmod(a, b, c) == 95367431640625u); 893 a = 100; b = 7919; c = 18446744073709551557u; 894 assert(powmod(a, b, c) == 18223853583554725198u); 895 a = 117; b = 7919; c = 18446744073709551557u; 896 assert(powmod(a, b, c) == 11493139548346411394u); 897 a = 134; b = 7919; c = 18446744073709551557u; 898 assert(powmod(a, b, c) == 10979163786734356774u); 899 a = 151; b = 7919; c = 18446744073709551557u; 900 assert(powmod(a, b, c) == 7023018419737782840u); 901 a = 168; b = 7919; c = 18446744073709551557u; 902 assert(powmod(a, b, c) == 58082701842386811u); 903 a = 185; b = 7919; c = 18446744073709551557u; 904 assert(powmod(a, b, c) == 17423478386299876798u); 905 a = 202; b = 7919; c = 18446744073709551557u; 906 assert(powmod(a, b, c) == 5522733478579799075u); 907 a = 219; b = 7919; c = 18446744073709551557u; 908 assert(powmod(a, b, c) == 15230218982491623487u); 909 a = 236; b = 7919; c = 18446744073709551557u; 910 assert(powmod(a, b, c) == 5198328724976436000u); 911 912 a = 0; b = 7919; c = 18446744073709551557u; 913 assert(powmod(a, b, c) == 0); 914 a = 123; b = 0; c = 18446744073709551557u; 915 assert(powmod(a, b, c) == 1); 916 917 immutable ulong a1 = 253, b1 = 7919, c1 = 18446744073709551557u; 918 assert(powmod(a1, b1, c1) == 3883707345459248860u); 919 920 uint x = 100 ,y = 7919, z = 1844674407u; 921 assert(powmod(x, y, z) == 1613100340u); 922 x = 134; y = 7919; z = 1844674407u; 923 assert(powmod(x, y, z) == 734956622u); 924 x = 151; y = 7919; z = 1844674407u; 925 assert(powmod(x, y, z) == 1738696945u); 926 x = 168; y = 7919; z = 1844674407u; 927 assert(powmod(x, y, z) == 1247580927u); 928 x = 185; y = 7919; z = 1844674407u; 929 assert(powmod(x, y, z) == 1293855176u); 930 x = 202; y = 7919; z = 1844674407u; 931 assert(powmod(x, y, z) == 1566963682u); 932 x = 219; y = 7919; z = 1844674407u; 933 assert(powmod(x, y, z) == 181227807u); 934 x = 236; y = 7919; z = 1844674407u; 935 assert(powmod(x, y, z) == 217988321u); 936 x = 253; y = 7919; z = 1844674407u; 937 assert(powmod(x, y, z) == 1588843243u); 938 939 x = 0; y = 7919; z = 184467u; 940 assert(powmod(x, y, z) == 0); 941 x = 123; y = 0; z = 1844674u; 942 assert(powmod(x, y, z) == 1); 943 944 immutable ubyte x1 = 117; 945 immutable uint y1 = 7919; 946 immutable uint z1 = 1844674407u; 947 auto res = powmod(x1, y1, z1); 948 assert(is(typeof(res) == uint)); 949 assert(res == 9479781u); 950 951 immutable ushort x2 = 123; 952 immutable uint y2 = 203; 953 immutable ubyte z2 = 113; 954 auto res2 = powmod(x2, y2, z2); 955 assert(is(typeof(res2) == ushort)); 956 assert(res2 == 42u); 957 } 958 959 /** 960 * Calculates e$(SUPERSCRIPT x). 961 * 962 * $(TABLE_SV 963 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)) ) 964 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) 965 * $(TR $(TD -$(INFIN)) $(TD +0.0) ) 966 * $(TR $(TD $(NAN)) $(TD $(NAN)) ) 967 * ) 968 */ 969 pragma(inline, true) 970 real exp(real x) @trusted pure nothrow @nogc // TODO: @safe 971 { 972 version (InlineAsm_X87) 973 { 974 import std.math.constants : LOG2E; 975 976 // e^^x = 2^^(LOG2E*x) 977 // (This is valid because the overflow & underflow limits for exp 978 // and exp2 are so similar). 979 if (!__ctfe) 980 return exp2Asm(LOG2E*x); 981 } 982 return expImpl(x); 983 } 984 985 /// ditto 986 pragma(inline, true) 987 double exp(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) exp(cast(real) x) : expImpl(x); } 988 989 /// ditto 990 pragma(inline, true) 991 float exp(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) exp(cast(real) x) : expImpl(x); } 992 993 /// 994 @safe unittest 995 { 996 import std.math.operations : feqrel; 997 import std.math.constants : E; 998 999 assert(exp(0.0) == 1.0); 1000 assert(exp(3.0).feqrel(E * E * E) > 16); 1001 } 1002 1003 private T expImpl(T)(T x) @safe pure nothrow @nogc 1004 { 1005 import std.math : floatTraits, RealFormat; 1006 import std.math.traits : isNaN; 1007 import std.math.rounding : floor; 1008 import std.math.algebraic : poly; 1009 import std.math.constants : LOG2E; 1010 1011 alias F = floatTraits!T; 1012 static if (F.realFormat == RealFormat.ieeeSingle) 1013 { 1014 static immutable T[6] P = [ 1015 5.0000001201E-1, 1016 1.6666665459E-1, 1017 4.1665795894E-2, 1018 8.3334519073E-3, 1019 1.3981999507E-3, 1020 1.9875691500E-4, 1021 ]; 1022 1023 enum T C1 = 0.693359375; 1024 enum T C2 = -2.12194440e-4; 1025 1026 // Overflow and Underflow limits. 1027 enum T OF = 88.72283905206835; 1028 enum T UF = -103.278929903431851103; // ln(2^-149) 1029 } 1030 else static if (F.realFormat == RealFormat.ieeeDouble) 1031 { 1032 // Coefficients for exp(x) 1033 static immutable T[3] P = [ 1034 9.99999999999999999910E-1L, 1035 3.02994407707441961300E-2L, 1036 1.26177193074810590878E-4L, 1037 ]; 1038 static immutable T[4] Q = [ 1039 2.00000000000000000009E0L, 1040 2.27265548208155028766E-1L, 1041 2.52448340349684104192E-3L, 1042 3.00198505138664455042E-6L, 1043 ]; 1044 1045 // C1 + C2 = LN2. 1046 enum T C1 = 6.93145751953125E-1; 1047 enum T C2 = 1.42860682030941723212E-6; 1048 1049 // Overflow and Underflow limits. 1050 enum T OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024) 1051 enum T UF = -7.451332191019412076235E2; // ln(2^-1075) 1052 } 1053 else static if (F.realFormat == RealFormat.ieeeExtended || 1054 F.realFormat == RealFormat.ieeeExtended53) 1055 { 1056 // Coefficients for exp(x) 1057 static immutable T[3] P = [ 1058 9.9999999999999999991025E-1L, 1059 3.0299440770744196129956E-2L, 1060 1.2617719307481059087798E-4L, 1061 ]; 1062 static immutable T[4] Q = [ 1063 2.0000000000000000000897E0L, 1064 2.2726554820815502876593E-1L, 1065 2.5244834034968410419224E-3L, 1066 3.0019850513866445504159E-6L, 1067 ]; 1068 1069 // C1 + C2 = LN2. 1070 enum T C1 = 6.9314575195312500000000E-1L; 1071 enum T C2 = 1.4286068203094172321215E-6L; 1072 1073 // Overflow and Underflow limits. 1074 enum T OF = 1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384) 1075 enum T UF = -1.13994985314888605586758E4L; // ln(2^-16446) 1076 } 1077 else static if (F.realFormat == RealFormat.ieeeQuadruple) 1078 { 1079 // Coefficients for exp(x) - 1 1080 static immutable T[5] P = [ 1081 9.999999999999999999999999999999999998502E-1L, 1082 3.508710990737834361215404761139478627390E-2L, 1083 2.708775201978218837374512615596512792224E-4L, 1084 6.141506007208645008909088812338454698548E-7L, 1085 3.279723985560247033712687707263393506266E-10L 1086 ]; 1087 static immutable T[6] Q = [ 1088 2.000000000000000000000000000000000000150E0, 1089 2.368408864814233538909747618894558968880E-1L, 1090 3.611828913847589925056132680618007270344E-3L, 1091 1.504792651814944826817779302637284053660E-5L, 1092 1.771372078166251484503904874657985291164E-8L, 1093 2.980756652081995192255342779918052538681E-12L 1094 ]; 1095 1096 // C1 + C2 = LN2. 1097 enum T C1 = 6.93145751953125E-1L; 1098 enum T C2 = 1.428606820309417232121458176568075500134E-6L; 1099 1100 // Overflow and Underflow limits. 1101 enum T OF = 1.135583025911358400418251384584930671458833e4L; 1102 enum T UF = -1.143276959615573793352782661133116431383730e4L; 1103 } 1104 else 1105 static assert(0, "Not implemented for this architecture"); 1106 1107 // Special cases. 1108 if (isNaN(x)) 1109 return x; 1110 if (x > OF) 1111 return real.infinity; 1112 if (x < UF) 1113 return 0.0; 1114 1115 // Express: e^^x = e^^g * 2^^n 1116 // = e^^g * e^^(n * LOG2E) 1117 // = e^^(g + n * LOG2E) 1118 T xx = floor((cast(T) LOG2E) * x + cast(T) 0.5); 1119 const int n = cast(int) xx; 1120 x -= xx * C1; 1121 x -= xx * C2; 1122 1123 static if (F.realFormat == RealFormat.ieeeSingle) 1124 { 1125 xx = x * x; 1126 x = poly(x, P) * xx + x + 1.0f; 1127 } 1128 else 1129 { 1130 // Rational approximation for exponential of the fractional part: 1131 // e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2)) 1132 xx = x * x; 1133 const T px = x * poly(xx, P); 1134 x = px / (poly(xx, Q) - px); 1135 x = (cast(T) 1.0) + (cast(T) 2.0) * x; 1136 } 1137 1138 // Scale by power of 2. 1139 x = core.math.ldexp(x, n); 1140 1141 return x; 1142 } 1143 1144 @safe @nogc nothrow unittest 1145 { 1146 import std.math : floatTraits, RealFormat; 1147 import std.math.operations : NaN, feqrel, isClose; 1148 import std.math.constants : E; 1149 import std.math.traits : isIdentical; 1150 import std.math.algebraic : abs; 1151 1152 version (IeeeFlagsSupport) import std.math.hardware : IeeeFlags, resetIeeeFlags, ieeeFlags; 1153 version (FloatingPointControlSupport) 1154 { 1155 import std.math.hardware : FloatingPointControl; 1156 1157 FloatingPointControl ctrl; 1158 if (FloatingPointControl.hasExceptionTraps) 1159 ctrl.disableExceptions(FloatingPointControl.allExceptions); 1160 ctrl.rounding = FloatingPointControl.roundToNearest; 1161 } 1162 1163 static void testExp(T)() 1164 { 1165 enum realFormat = floatTraits!T.realFormat; 1166 static if (realFormat == RealFormat.ieeeQuadruple) 1167 { 1168 static immutable T[2][] exptestpoints = 1169 [ // x exp(x) 1170 [ 1.0L, E ], 1171 [ 0.5L, 0x1.a61298e1e069bc972dfefab6df34p+0L ], 1172 [ 3.0L, E*E*E ], 1173 [ 0x1.6p+13L, 0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow 1174 [ 0x1.7p+13L, T.infinity ], // close overflow 1175 [ 0x1p+80L, T.infinity ], // far overflow 1176 [ T.infinity, T.infinity ], 1177 [-0x1.18p+13L, 0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow 1178 [-0x1.625p+13L, 0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto 1179 [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal 1180 [-0x1.6549p+13L, 0x0.0000000000000000000000000001p-16382L ], // ditto 1181 [-0x1.655p+13L, 0 ], // close underflow 1182 [-0x1p+30L, 0 ], // far underflow 1183 ]; 1184 } 1185 else static if (realFormat == RealFormat.ieeeExtended || 1186 realFormat == RealFormat.ieeeExtended53) 1187 { 1188 static immutable T[2][] exptestpoints = 1189 [ // x exp(x) 1190 [ 1.0L, E ], 1191 [ 0.5L, 0x1.a61298e1e069bc97p+0L ], 1192 [ 3.0L, E*E*E ], 1193 [ 0x1.1p+13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow 1194 [ 0x1.7p+13L, T.infinity ], // close overflow 1195 [ 0x1p+80L, T.infinity ], // far overflow 1196 [ T.infinity, T.infinity ], 1197 [-0x1.18p+13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow 1198 [-0x1.625p+13L, 0x1.a6bd68a39d11f35cp-16358L ], // ditto 1199 [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L ], // near underflow - subnormal 1200 [-0x1.643p+13L, 0x1p-16444L ], // ditto 1201 [-0x1.645p+13L, 0 ], // close underflow 1202 [-0x1p+30L, 0 ], // far underflow 1203 ]; 1204 } 1205 else static if (realFormat == RealFormat.ieeeDouble) 1206 { 1207 static immutable T[2][] exptestpoints = 1208 [ // x, exp(x) 1209 [ 1.0L, E ], 1210 [ 0.5L, 0x1.a61298e1e069cp+0L ], 1211 [ 3.0L, E*E*E ], 1212 [ 0x1.6p+9L, 0x1.93bf4ec282efbp+1015L ], // near overflow 1213 [ 0x1.7p+9L, T.infinity ], // close overflow 1214 [ 0x1p+80L, T.infinity ], // far overflow 1215 [ T.infinity, T.infinity ], 1216 [-0x1.6p+9L, 0x1.44a3824e5285fp-1016L ], // near underflow 1217 [-0x1.64p+9L, 0x0.06f84920bb2d4p-1022L ], // near underflow - subnormal 1218 [-0x1.743p+9L, 0x0.0000000000001p-1022L ], // ditto 1219 [-0x1.8p+9L, 0 ], // close underflow 1220 [-0x1p+30L, 0 ], // far underflow 1221 ]; 1222 } 1223 else static if (realFormat == RealFormat.ieeeSingle) 1224 { 1225 static immutable T[2][] exptestpoints = 1226 [ // x, exp(x) 1227 [ 1.0L, E ], 1228 [ 0.5L, 0x1.a61299p+0L ], 1229 [ 3.0L, E*E*E ], 1230 [ 0x1.62p+6L, 0x1.99b988p+127L ], // near overflow 1231 [ 0x1.7p+6L, T.infinity ], // close overflow 1232 [ 0x1p+80L, T.infinity ], // far overflow 1233 [ T.infinity, T.infinity ], 1234 [-0x1.5cp+6L, 0x1.666d0ep-126L ], // near underflow 1235 [-0x1.7p+6L, 0x0.026a42p-126L ], // near underflow - subnormal 1236 [-0x1.9cp+6L, 0x0.000002p-126L ], // ditto 1237 [-0x1.ap+6L, 0 ], // close underflow 1238 [-0x1p+30L, 0 ], // far underflow 1239 ]; 1240 } 1241 else 1242 static assert(0, "No exp() tests for real type!"); 1243 1244 const minEqualMantissaBits = T.mant_dig - 2; 1245 T x; 1246 version (IeeeFlagsSupport) IeeeFlags f; 1247 foreach (ref pair; exptestpoints) 1248 { 1249 version (IeeeFlagsSupport) resetIeeeFlags(); 1250 x = exp(pair[0]); 1251 //printf("exp(%La) = %La, should be %La\n", cast(real) pair[0], cast(real) x, cast(real) pair[1]); 1252 assert(feqrel(x, pair[1]) >= minEqualMantissaBits); 1253 } 1254 1255 // Ideally, exp(0) would not set the inexact flag. 1256 // Unfortunately, fldl2e sets it! 1257 // So it's not realistic to avoid setting it. 1258 assert(exp(cast(T) 0.0) == 1.0); 1259 1260 // NaN propagation. Doesn't set flags, bcos was already NaN. 1261 version (IeeeFlagsSupport) 1262 { 1263 resetIeeeFlags(); 1264 x = exp(T.nan); 1265 f = ieeeFlags; 1266 assert(isIdentical(abs(x), T.nan)); 1267 assert(f.flags == 0); 1268 1269 resetIeeeFlags(); 1270 x = exp(-T.nan); 1271 f = ieeeFlags; 1272 assert(isIdentical(abs(x), T.nan)); 1273 assert(f.flags == 0); 1274 } 1275 else 1276 { 1277 x = exp(T.nan); 1278 assert(isIdentical(abs(x), T.nan)); 1279 1280 x = exp(-T.nan); 1281 assert(isIdentical(abs(x), T.nan)); 1282 } 1283 1284 x = exp(NaN(0x123)); 1285 assert(isIdentical(x, NaN(0x123))); 1286 } 1287 1288 import std.meta : AliasSeq; 1289 foreach (T; AliasSeq!(real, double, float)) 1290 testExp!T(); 1291 1292 // High resolution test (verified against GNU MPFR/Mathematica). 1293 assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+0L); 1294 1295 assert(isClose(exp(3.0L), E * E * E, real.sizeof > double.sizeof ? 1e-15 : 1e-14)); 1296 } 1297 1298 /** 1299 * Calculates the value of the natural logarithm base (e) 1300 * raised to the power of x, minus 1. 1301 * 1302 * For very small x, expm1(x) is more accurate 1303 * than exp(x)-1. 1304 * 1305 * $(TABLE_SV 1306 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)-1) ) 1307 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) ) 1308 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) 1309 * $(TR $(TD -$(INFIN)) $(TD -1.0) ) 1310 * $(TR $(TD $(NAN)) $(TD $(NAN)) ) 1311 * ) 1312 */ 1313 pragma(inline, true) 1314 real expm1(real x) @trusted pure nothrow @nogc // TODO: @safe 1315 { 1316 version (InlineAsm_X87) 1317 { 1318 if (!__ctfe) 1319 return expm1Asm(x); 1320 } 1321 return expm1Impl(x); 1322 } 1323 1324 /// ditto 1325 pragma(inline, true) 1326 double expm1(double x) @safe pure nothrow @nogc 1327 { 1328 return __ctfe ? cast(double) expm1(cast(real) x) : expm1Impl(x); 1329 } 1330 1331 /// ditto 1332 pragma(inline, true) 1333 float expm1(float x) @safe pure nothrow @nogc 1334 { 1335 // no single-precision version in Cephes => use double precision 1336 return __ctfe ? cast(float) expm1(cast(real) x) : cast(float) expm1Impl(cast(double) x); 1337 } 1338 1339 /// 1340 @safe unittest 1341 { 1342 import std.math.traits : isIdentical; 1343 import std.math.operations : feqrel; 1344 1345 assert(isIdentical(expm1(0.0), 0.0)); 1346 assert(expm1(1.0).feqrel(1.71828) > 16); 1347 assert(expm1(2.0).feqrel(6.3890) > 16); 1348 } 1349 1350 version (InlineAsm_X87) 1351 private real expm1Asm(real x) @trusted pure nothrow @nogc 1352 { 1353 version (X86) 1354 { 1355 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4 1356 asm pure nothrow @nogc 1357 { 1358 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. 1359 * Author: Don Clugston. 1360 * 1361 * expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x. 1362 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y)) 1363 * and 2ym1 = (2^^(y-rndint(y))-1). 1364 * If 2rndy < 0.5*real.epsilon, result is -1. 1365 * Implementation is otherwise the same as for exp2() 1366 */ 1367 naked; 1368 fld real ptr [ESP+4] ; // x 1369 mov AX, [ESP+4+8]; // AX = exponent and sign 1370 sub ESP, 12+8; // Create scratch space on the stack 1371 // [ESP,ESP+2] = scratchint 1372 // [ESP+4..+6, +8..+10, +10] = scratchreal 1373 // set scratchreal mantissa = 1.0 1374 mov dword ptr [ESP+8], 0; 1375 mov dword ptr [ESP+8+4], 0x80000000; 1376 and AX, 0x7FFF; // drop sign bit 1377 cmp AX, 0x401D; // avoid InvalidException in fist 1378 jae L_extreme; 1379 fldl2e; 1380 fmulp ST(1), ST; // y = x*log2(e) 1381 fist dword ptr [ESP]; // scratchint = rndint(y) 1382 fisub dword ptr [ESP]; // y - rndint(y) 1383 // and now set scratchreal exponent 1384 mov EAX, [ESP]; 1385 add EAX, 0x3fff; 1386 jle short L_largenegative; 1387 cmp EAX,0x8000; 1388 jge short L_largepositive; 1389 mov [ESP+8+8],AX; 1390 f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1 1391 fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y) 1392 fmul ST(1), ST; // ST=2rndy, ST(1)=2rndy*2ym1 1393 fld1; 1394 fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1 1395 faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1 1396 add ESP,12+8; 1397 ret PARAMSIZE; 1398 1399 L_extreme: // Extreme exponent. X is very large positive, very 1400 // large negative, infinity, or NaN. 1401 fxam; 1402 fstsw AX; 1403 test AX, 0x0400; // NaN_or_zero, but we already know x != 0 1404 jz L_was_nan; // if x is NaN, returns x 1405 test AX, 0x0200; 1406 jnz L_largenegative; 1407 L_largepositive: 1408 // Set scratchreal = real.max. 1409 // squaring it will create infinity, and set overflow flag. 1410 mov word ptr [ESP+8+8], 0x7FFE; 1411 fstp ST(0); 1412 fld real ptr [ESP+8]; // load scratchreal 1413 fmul ST(0), ST; // square it, to create havoc! 1414 L_was_nan: 1415 add ESP,12+8; 1416 ret PARAMSIZE; 1417 L_largenegative: 1418 fstp ST(0); 1419 fld1; 1420 fchs; // return -1. Underflow flag is not set. 1421 add ESP,12+8; 1422 ret PARAMSIZE; 1423 } 1424 } 1425 else version (X86_64) 1426 { 1427 asm pure nothrow @nogc 1428 { 1429 naked; 1430 } 1431 version (Win64) 1432 { 1433 asm pure nothrow @nogc 1434 { 1435 fld real ptr [RCX]; // x 1436 mov AX,[RCX+8]; // AX = exponent and sign 1437 } 1438 } 1439 else 1440 { 1441 asm pure nothrow @nogc 1442 { 1443 fld real ptr [RSP+8]; // x 1444 mov AX,[RSP+8+8]; // AX = exponent and sign 1445 } 1446 } 1447 asm pure nothrow @nogc 1448 { 1449 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. 1450 * Author: Don Clugston. 1451 * 1452 * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x. 1453 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y)) 1454 * and 2ym1 = (2^(y-rndint(y))-1). 1455 * If 2rndy < 0.5*real.epsilon, result is -1. 1456 * Implementation is otherwise the same as for exp2() 1457 */ 1458 sub RSP, 24; // Create scratch space on the stack 1459 // [RSP,RSP+2] = scratchint 1460 // [RSP+4..+6, +8..+10, +10] = scratchreal 1461 // set scratchreal mantissa = 1.0 1462 mov dword ptr [RSP+8], 0; 1463 mov dword ptr [RSP+8+4], 0x80000000; 1464 and AX, 0x7FFF; // drop sign bit 1465 cmp AX, 0x401D; // avoid InvalidException in fist 1466 jae L_extreme; 1467 fldl2e; 1468 fmul ; // y = x*log2(e) 1469 fist dword ptr [RSP]; // scratchint = rndint(y) 1470 fisub dword ptr [RSP]; // y - rndint(y) 1471 // and now set scratchreal exponent 1472 mov EAX, [RSP]; 1473 add EAX, 0x3fff; 1474 jle short L_largenegative; 1475 cmp EAX,0x8000; 1476 jge short L_largepositive; 1477 mov [RSP+8+8],AX; 1478 f2xm1; // 2^(y-rndint(y)) -1 1479 fld real ptr [RSP+8] ; // 2^rndint(y) 1480 fmul ST(1), ST; 1481 fld1; 1482 fsubp ST(1), ST; 1483 fadd; 1484 add RSP,24; 1485 ret; 1486 1487 L_extreme: // Extreme exponent. X is very large positive, very 1488 // large negative, infinity, or NaN. 1489 fxam; 1490 fstsw AX; 1491 test AX, 0x0400; // NaN_or_zero, but we already know x != 0 1492 jz L_was_nan; // if x is NaN, returns x 1493 test AX, 0x0200; 1494 jnz L_largenegative; 1495 L_largepositive: 1496 // Set scratchreal = real.max. 1497 // squaring it will create infinity, and set overflow flag. 1498 mov word ptr [RSP+8+8], 0x7FFE; 1499 fstp ST(0); 1500 fld real ptr [RSP+8]; // load scratchreal 1501 fmul ST(0), ST; // square it, to create havoc! 1502 L_was_nan: 1503 add RSP,24; 1504 ret; 1505 1506 L_largenegative: 1507 fstp ST(0); 1508 fld1; 1509 fchs; // return -1. Underflow flag is not set. 1510 add RSP,24; 1511 ret; 1512 } 1513 } 1514 else 1515 static assert(0); 1516 } 1517 1518 private T expm1Impl(T)(T x) @safe pure nothrow @nogc 1519 { 1520 import std.math : floatTraits, RealFormat; 1521 import std.math.rounding : floor; 1522 import std.math.algebraic : poly; 1523 import std.math.constants : LN2; 1524 1525 // Coefficients for exp(x) - 1 and overflow/underflow limits. 1526 enum realFormat = floatTraits!T.realFormat; 1527 static if (realFormat == RealFormat.ieeeQuadruple) 1528 { 1529 static immutable T[8] P = [ 1530 2.943520915569954073888921213330863757240E8L, 1531 -5.722847283900608941516165725053359168840E7L, 1532 8.944630806357575461578107295909719817253E6L, 1533 -7.212432713558031519943281748462837065308E5L, 1534 4.578962475841642634225390068461943438441E4L, 1535 -1.716772506388927649032068540558788106762E3L, 1536 4.401308817383362136048032038528753151144E1L, 1537 -4.888737542888633647784737721812546636240E-1L 1538 ]; 1539 1540 static immutable T[9] Q = [ 1541 1.766112549341972444333352727998584753865E9L, 1542 -7.848989743695296475743081255027098295771E8L, 1543 1.615869009634292424463780387327037251069E8L, 1544 -2.019684072836541751428967854947019415698E7L, 1545 1.682912729190313538934190635536631941751E6L, 1546 -9.615511549171441430850103489315371768998E4L, 1547 3.697714952261803935521187272204485251835E3L, 1548 -8.802340681794263968892934703309274564037E1L, 1549 1.0 1550 ]; 1551 1552 enum T OF = 1.1356523406294143949491931077970764891253E4L; 1553 enum T UF = -1.143276959615573793352782661133116431383730e4L; 1554 } 1555 else static if (realFormat == RealFormat.ieeeExtended) 1556 { 1557 static immutable T[5] P = [ 1558 -1.586135578666346600772998894928250240826E4L, 1559 2.642771505685952966904660652518429479531E3L, 1560 -3.423199068835684263987132888286791620673E2L, 1561 1.800826371455042224581246202420972737840E1L, 1562 -5.238523121205561042771939008061958820811E-1L, 1563 ]; 1564 static immutable T[6] Q = [ 1565 -9.516813471998079611319047060563358064497E4L, 1566 3.964866271411091674556850458227710004570E4L, 1567 -7.207678383830091850230366618190187434796E3L, 1568 7.206038318724600171970199625081491823079E2L, 1569 -4.002027679107076077238836622982900945173E1L, 1570 1.0 1571 ]; 1572 1573 enum T OF = 1.1356523406294143949492E4L; 1574 enum T UF = -4.5054566736396445112120088E1L; 1575 } 1576 else static if (realFormat == RealFormat.ieeeDouble) 1577 { 1578 static immutable T[3] P = [ 1579 9.9999999999999999991025E-1, 1580 3.0299440770744196129956E-2, 1581 1.2617719307481059087798E-4, 1582 ]; 1583 static immutable T[4] Q = [ 1584 2.0000000000000000000897E0, 1585 2.2726554820815502876593E-1, 1586 2.5244834034968410419224E-3, 1587 3.0019850513866445504159E-6, 1588 ]; 1589 } 1590 else 1591 static assert(0, "no coefficients for expm1()"); 1592 1593 static if (realFormat == RealFormat.ieeeDouble) // special case for double precision 1594 { 1595 if (x < -0.5 || x > 0.5) 1596 return exp(x) - 1.0; 1597 if (x == 0.0) 1598 return x; 1599 1600 const T xx = x * x; 1601 x = x * poly(xx, P); 1602 x = x / (poly(xx, Q) - x); 1603 return x + x; 1604 } 1605 else 1606 { 1607 // C1 + C2 = LN2. 1608 enum T C1 = 6.9314575195312500000000E-1L; 1609 enum T C2 = 1.428606820309417232121458176568075500134E-6L; 1610 1611 // Special cases. 1612 if (x > OF) 1613 return real.infinity; 1614 if (x == cast(T) 0.0) 1615 return x; 1616 if (x < UF) 1617 return -1.0; 1618 1619 // Express x = LN2 (n + remainder), remainder not exceeding 1/2. 1620 int n = cast(int) floor((cast(T) 0.5) + x / cast(T) LN2); 1621 x -= n * C1; 1622 x -= n * C2; 1623 1624 // Rational approximation: 1625 // exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x) 1626 T px = x * poly(x, P); 1627 T qx = poly(x, Q); 1628 const T xx = x * x; 1629 qx = x + ((cast(T) 0.5) * xx + xx * px / qx); 1630 1631 // We have qx = exp(remainder LN2) - 1, so: 1632 // exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1. 1633 px = core.math.ldexp(cast(T) 1.0, n); 1634 x = px * qx + (px - cast(T) 1.0); 1635 1636 return x; 1637 } 1638 } 1639 1640 @safe @nogc nothrow unittest 1641 { 1642 import std.math.traits : isNaN; 1643 import std.math.operations : isClose, CommonDefaultFor; 1644 1645 static void testExpm1(T)() 1646 { 1647 // NaN 1648 assert(isNaN(expm1(cast(T) T.nan))); 1649 1650 static immutable T[] xs = [ -2, -0.75, -0.3, 0.0, 0.1, 0.2, 0.5, 1.0 ]; 1651 foreach (x; xs) 1652 { 1653 const T e = expm1(x); 1654 const T r = exp(x) - 1; 1655 1656 //printf("expm1(%Lg) = %Lg, should approximately be %Lg\n", cast(real) x, cast(real) e, cast(real) r); 1657 assert(isClose(r, e, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T))); 1658 } 1659 } 1660 1661 import std.meta : AliasSeq; 1662 foreach (T; AliasSeq!(real, double)) 1663 testExpm1!T(); 1664 } 1665 1666 /** 1667 * Calculates 2$(SUPERSCRIPT x). 1668 * 1669 * $(TABLE_SV 1670 * $(TR $(TH x) $(TH exp2(x)) ) 1671 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) 1672 * $(TR $(TD -$(INFIN)) $(TD +0.0) ) 1673 * $(TR $(TD $(NAN)) $(TD $(NAN)) ) 1674 * ) 1675 */ 1676 pragma(inline, true) 1677 real exp2(real x) @nogc @trusted pure nothrow // TODO: @safe 1678 { 1679 version (InlineAsm_X87) 1680 { 1681 if (!__ctfe) 1682 return exp2Asm(x); 1683 } 1684 return exp2Impl(x); 1685 } 1686 1687 /// ditto 1688 pragma(inline, true) 1689 double exp2(double x) @nogc @safe pure nothrow { return __ctfe ? cast(double) exp2(cast(real) x) : exp2Impl(x); } 1690 1691 /// ditto 1692 pragma(inline, true) 1693 float exp2(float x) @nogc @safe pure nothrow { return __ctfe ? cast(float) exp2(cast(real) x) : exp2Impl(x); } 1694 1695 /// 1696 @safe unittest 1697 { 1698 import std.math.traits : isIdentical; 1699 import std.math.operations : feqrel; 1700 1701 assert(isIdentical(exp2(0.0), 1.0)); 1702 assert(exp2(2.0).feqrel(4.0) > 16); 1703 assert(exp2(8.0).feqrel(256.0) > 16); 1704 } 1705 1706 @safe unittest 1707 { 1708 version (CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented 1709 { 1710 assert( core.stdc.math.exp2f(0.0f) == 1 ); 1711 assert( core.stdc.math.exp2 (0.0) == 1 ); 1712 assert( core.stdc.math.exp2l(0.0L) == 1 ); 1713 } 1714 } 1715 1716 version (InlineAsm_X87) 1717 private real exp2Asm(real x) @nogc @trusted pure nothrow 1718 { 1719 version (X86) 1720 { 1721 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4 1722 1723 asm pure nothrow @nogc 1724 { 1725 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. 1726 * Author: Don Clugston. 1727 * 1728 * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x)) 1729 * The trick for high performance is to avoid the fscale(28cycles on core2), 1730 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. 1731 * 1732 * We can do frndint by using fist. BUT we can't use it for huge numbers, 1733 * because it will set the Invalid Operation flag if overflow or NaN occurs. 1734 * Fortunately, whenever this happens the result would be zero or infinity. 1735 * 1736 * We can perform fscale by directly poking into the exponent. BUT this doesn't 1737 * work for the (very rare) cases where the result is subnormal. So we fall back 1738 * to the slow method in that case. 1739 */ 1740 naked; 1741 fld real ptr [ESP+4] ; // x 1742 mov AX, [ESP+4+8]; // AX = exponent and sign 1743 sub ESP, 12+8; // Create scratch space on the stack 1744 // [ESP,ESP+2] = scratchint 1745 // [ESP+4..+6, +8..+10, +10] = scratchreal 1746 // set scratchreal mantissa = 1.0 1747 mov dword ptr [ESP+8], 0; 1748 mov dword ptr [ESP+8+4], 0x80000000; 1749 and AX, 0x7FFF; // drop sign bit 1750 cmp AX, 0x401D; // avoid InvalidException in fist 1751 jae L_extreme; 1752 fist dword ptr [ESP]; // scratchint = rndint(x) 1753 fisub dword ptr [ESP]; // x - rndint(x) 1754 // and now set scratchreal exponent 1755 mov EAX, [ESP]; 1756 add EAX, 0x3fff; 1757 jle short L_subnormal; 1758 cmp EAX,0x8000; 1759 jge short L_overflow; 1760 mov [ESP+8+8],AX; 1761 L_normal: 1762 f2xm1; 1763 fld1; 1764 faddp ST(1), ST; // 2^^(x-rndint(x)) 1765 fld real ptr [ESP+8] ; // 2^^rndint(x) 1766 add ESP,12+8; 1767 fmulp ST(1), ST; 1768 ret PARAMSIZE; 1769 1770 L_subnormal: 1771 // Result will be subnormal. 1772 // In this rare case, the simple poking method doesn't work. 1773 // The speed doesn't matter, so use the slow fscale method. 1774 fild dword ptr [ESP]; // scratchint 1775 fld1; 1776 fscale; 1777 fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint 1778 fstp ST(0); // drop scratchint 1779 jmp L_normal; 1780 1781 L_extreme: // Extreme exponent. X is very large positive, very 1782 // large negative, infinity, or NaN. 1783 fxam; 1784 fstsw AX; 1785 test AX, 0x0400; // NaN_or_zero, but we already know x != 0 1786 jz L_was_nan; // if x is NaN, returns x 1787 // set scratchreal = real.min_normal 1788 // squaring it will return 0, setting underflow flag 1789 mov word ptr [ESP+8+8], 1; 1790 test AX, 0x0200; 1791 jnz L_waslargenegative; 1792 L_overflow: 1793 // Set scratchreal = real.max. 1794 // squaring it will create infinity, and set overflow flag. 1795 mov word ptr [ESP+8+8], 0x7FFE; 1796 L_waslargenegative: 1797 fstp ST(0); 1798 fld real ptr [ESP+8]; // load scratchreal 1799 fmul ST(0), ST; // square it, to create havoc! 1800 L_was_nan: 1801 add ESP,12+8; 1802 ret PARAMSIZE; 1803 } 1804 } 1805 else version (X86_64) 1806 { 1807 asm pure nothrow @nogc 1808 { 1809 naked; 1810 } 1811 version (Win64) 1812 { 1813 asm pure nothrow @nogc 1814 { 1815 fld real ptr [RCX]; // x 1816 mov AX,[RCX+8]; // AX = exponent and sign 1817 } 1818 } 1819 else 1820 { 1821 asm pure nothrow @nogc 1822 { 1823 fld real ptr [RSP+8]; // x 1824 mov AX,[RSP+8+8]; // AX = exponent and sign 1825 } 1826 } 1827 asm pure nothrow @nogc 1828 { 1829 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. 1830 * Author: Don Clugston. 1831 * 1832 * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x)) 1833 * The trick for high performance is to avoid the fscale(28cycles on core2), 1834 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. 1835 * 1836 * We can do frndint by using fist. BUT we can't use it for huge numbers, 1837 * because it will set the Invalid Operation flag is overflow or NaN occurs. 1838 * Fortunately, whenever this happens the result would be zero or infinity. 1839 * 1840 * We can perform fscale by directly poking into the exponent. BUT this doesn't 1841 * work for the (very rare) cases where the result is subnormal. So we fall back 1842 * to the slow method in that case. 1843 */ 1844 sub RSP, 24; // Create scratch space on the stack 1845 // [RSP,RSP+2] = scratchint 1846 // [RSP+4..+6, +8..+10, +10] = scratchreal 1847 // set scratchreal mantissa = 1.0 1848 mov dword ptr [RSP+8], 0; 1849 mov dword ptr [RSP+8+4], 0x80000000; 1850 and AX, 0x7FFF; // drop sign bit 1851 cmp AX, 0x401D; // avoid InvalidException in fist 1852 jae L_extreme; 1853 fist dword ptr [RSP]; // scratchint = rndint(x) 1854 fisub dword ptr [RSP]; // x - rndint(x) 1855 // and now set scratchreal exponent 1856 mov EAX, [RSP]; 1857 add EAX, 0x3fff; 1858 jle short L_subnormal; 1859 cmp EAX,0x8000; 1860 jge short L_overflow; 1861 mov [RSP+8+8],AX; 1862 L_normal: 1863 f2xm1; 1864 fld1; 1865 fadd; // 2^(x-rndint(x)) 1866 fld real ptr [RSP+8] ; // 2^rndint(x) 1867 add RSP,24; 1868 fmulp ST(1), ST; 1869 ret; 1870 1871 L_subnormal: 1872 // Result will be subnormal. 1873 // In this rare case, the simple poking method doesn't work. 1874 // The speed doesn't matter, so use the slow fscale method. 1875 fild dword ptr [RSP]; // scratchint 1876 fld1; 1877 fscale; 1878 fstp real ptr [RSP+8]; // scratchreal = 2^scratchint 1879 fstp ST(0); // drop scratchint 1880 jmp L_normal; 1881 1882 L_extreme: // Extreme exponent. X is very large positive, very 1883 // large negative, infinity, or NaN. 1884 fxam; 1885 fstsw AX; 1886 test AX, 0x0400; // NaN_or_zero, but we already know x != 0 1887 jz L_was_nan; // if x is NaN, returns x 1888 // set scratchreal = real.min 1889 // squaring it will return 0, setting underflow flag 1890 mov word ptr [RSP+8+8], 1; 1891 test AX, 0x0200; 1892 jnz L_waslargenegative; 1893 L_overflow: 1894 // Set scratchreal = real.max. 1895 // squaring it will create infinity, and set overflow flag. 1896 mov word ptr [RSP+8+8], 0x7FFE; 1897 L_waslargenegative: 1898 fstp ST(0); 1899 fld real ptr [RSP+8]; // load scratchreal 1900 fmul ST(0), ST; // square it, to create havoc! 1901 L_was_nan: 1902 add RSP,24; 1903 ret; 1904 } 1905 } 1906 else 1907 static assert(0); 1908 } 1909 1910 private T exp2Impl(T)(T x) @nogc @safe pure nothrow 1911 { 1912 import std.math : floatTraits, RealFormat; 1913 import std.math.traits : isNaN; 1914 import std.math.rounding : floor; 1915 import std.math.algebraic : poly; 1916 1917 // Coefficients for exp2(x) 1918 enum realFormat = floatTraits!T.realFormat; 1919 static if (realFormat == RealFormat.ieeeQuadruple) 1920 { 1921 static immutable T[5] P = [ 1922 9.079594442980146270952372234833529694788E12L, 1923 1.530625323728429161131811299626419117557E11L, 1924 5.677513871931844661829755443994214173883E8L, 1925 6.185032670011643762127954396427045467506E5L, 1926 1.587171580015525194694938306936721666031E2L 1927 ]; 1928 1929 static immutable T[6] Q = [ 1930 2.619817175234089411411070339065679229869E13L, 1931 1.490560994263653042761789432690793026977E12L, 1932 1.092141473886177435056423606755843616331E10L, 1933 2.186249607051644894762167991800811827835E7L, 1934 1.236602014442099053716561665053645270207E4L, 1935 1.0 1936 ]; 1937 } 1938 else static if (realFormat == RealFormat.ieeeExtended) 1939 { 1940 static immutable T[3] P = [ 1941 2.0803843631901852422887E6L, 1942 3.0286971917562792508623E4L, 1943 6.0614853552242266094567E1L, 1944 ]; 1945 static immutable T[4] Q = [ 1946 6.0027204078348487957118E6L, 1947 3.2772515434906797273099E5L, 1948 1.7492876999891839021063E3L, 1949 1.0000000000000000000000E0L, 1950 ]; 1951 } 1952 else static if (realFormat == RealFormat.ieeeDouble) 1953 { 1954 static immutable T[3] P = [ 1955 1.51390680115615096133E3L, 1956 2.02020656693165307700E1L, 1957 2.30933477057345225087E-2L, 1958 ]; 1959 static immutable T[3] Q = [ 1960 4.36821166879210612817E3L, 1961 2.33184211722314911771E2L, 1962 1.00000000000000000000E0L, 1963 ]; 1964 } 1965 else static if (realFormat == RealFormat.ieeeSingle) 1966 { 1967 static immutable T[6] P = [ 1968 6.931472028550421E-001L, 1969 2.402264791363012E-001L, 1970 5.550332471162809E-002L, 1971 9.618437357674640E-003L, 1972 1.339887440266574E-003L, 1973 1.535336188319500E-004L, 1974 ]; 1975 } 1976 else 1977 static assert(0, "no coefficients for exp2()"); 1978 1979 // Overflow and Underflow limits. 1980 enum T OF = T.max_exp; 1981 enum T UF = T.min_exp - 1; 1982 1983 // Special cases. 1984 if (isNaN(x)) 1985 return x; 1986 if (x > OF) 1987 return real.infinity; 1988 if (x < UF) 1989 return 0.0; 1990 1991 static if (realFormat == RealFormat.ieeeSingle) // special case for single precision 1992 { 1993 // The following is necessary because range reduction blows up. 1994 if (x == 0.0f) 1995 return 1.0f; 1996 1997 // Separate into integer and fractional parts. 1998 const T i = floor(x); 1999 int n = cast(int) i; 2000 x -= i; 2001 if (x > 0.5f) 2002 { 2003 n += 1; 2004 x -= 1.0f; 2005 } 2006 2007 // Rational approximation: 2008 // exp2(x) = 1.0 + x P(x) 2009 x = 1.0f + x * poly(x, P); 2010 } 2011 else 2012 { 2013 // Separate into integer and fractional parts. 2014 const T i = floor(x + cast(T) 0.5); 2015 int n = cast(int) i; 2016 x -= i; 2017 2018 // Rational approximation: 2019 // exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2)) 2020 const T xx = x * x; 2021 const T px = x * poly(xx, P); 2022 x = px / (poly(xx, Q) - px); 2023 x = (cast(T) 1.0) + (cast(T) 2.0) * x; 2024 } 2025 2026 // Scale by power of 2. 2027 x = core.math.ldexp(x, n); 2028 2029 return x; 2030 } 2031 2032 @safe @nogc nothrow unittest 2033 { 2034 import std.math.operations : feqrel, NaN, isClose; 2035 import std.math.traits : isIdentical; 2036 import std.math.constants : SQRT2; 2037 2038 assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1); 2039 assert(exp2(8.0L) == 256.0); 2040 assert(exp2(-9.0L)== 1.0L/512.0); 2041 2042 static void testExp2(T)() 2043 { 2044 // NaN 2045 const T specialNaN = NaN(0x0123L); 2046 assert(isIdentical(exp2(specialNaN), specialNaN)); 2047 2048 // over-/underflow 2049 enum T OF = T.max_exp; 2050 enum T UF = T.min_exp - T.mant_dig; 2051 assert(isIdentical(exp2(OF + 1), cast(T) T.infinity)); 2052 assert(isIdentical(exp2(UF - 1), cast(T) 0.0)); 2053 2054 static immutable T[2][] vals = 2055 [ 2056 // x, exp2(x) 2057 [ 0.0, 1.0 ], 2058 [ -0.0, 1.0 ], 2059 [ 0.5, SQRT2 ], 2060 [ 8.0, 256.0 ], 2061 [ -9.0, 1.0 / 512 ], 2062 ]; 2063 2064 foreach (ref val; vals) 2065 { 2066 const T x = val[0]; 2067 const T r = val[1]; 2068 const T e = exp2(x); 2069 2070 //printf("exp2(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) e, cast(real) r); 2071 assert(isClose(r, e)); 2072 } 2073 } 2074 2075 import std.meta : AliasSeq; 2076 foreach (T; AliasSeq!(real, double, float)) 2077 testExp2!T(); 2078 } 2079 2080 /********************************************************************* 2081 * Separate floating point value into significand and exponent. 2082 * 2083 * Returns: 2084 * Calculate and return $(I x) and $(I exp) such that 2085 * value =$(I x)*2$(SUPERSCRIPT exp) and 2086 * .5 $(LT)= |$(I x)| $(LT) 1.0 2087 * 2088 * $(I x) has same sign as value. 2089 * 2090 * $(TABLE_SV 2091 * $(TR $(TH value) $(TH returns) $(TH exp)) 2092 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0)) 2093 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max)) 2094 * $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min)) 2095 * $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min)) 2096 * ) 2097 */ 2098 T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc 2099 if (isFloatingPoint!T) 2100 { 2101 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB; 2102 import std.math.traits : isSubnormal; 2103 2104 if (__ctfe) 2105 { 2106 // Handle special cases. 2107 if (value == 0) { exp = 0; return value; } 2108 if (value == T.infinity) { exp = int.max; return value; } 2109 if (value == -T.infinity || value != value) { exp = int.min; return value; } 2110 // Handle ordinary cases. 2111 // In CTFE there is no performance advantage for having separate 2112 // paths for different floating point types. 2113 T absValue = value < 0 ? -value : value; 2114 int expCount; 2115 static if (T.mant_dig > double.mant_dig) 2116 { 2117 for (; absValue >= 0x1.0p+1024L; absValue *= 0x1.0p-1024L) 2118 expCount += 1024; 2119 for (; absValue < 0x1.0p-1021L; absValue *= 0x1.0p+1021L) 2120 expCount -= 1021; 2121 } 2122 const double dval = cast(double) absValue; 2123 int dexp = cast(int) (((*cast(const long*) &dval) >>> 52) & 0x7FF) + double.min_exp - 2; 2124 dexp++; 2125 expCount += dexp; 2126 absValue *= 2.0 ^^ -dexp; 2127 // If the original value was subnormal or if it was a real 2128 // then absValue can still be outside the [0.5, 1.0) range. 2129 if (absValue < 0.5) 2130 { 2131 assert(T.mant_dig > double.mant_dig || isSubnormal(value)); 2132 do 2133 { 2134 absValue += absValue; 2135 expCount--; 2136 } while (absValue < 0.5); 2137 } 2138 else 2139 { 2140 assert(absValue < 1 || T.mant_dig > double.mant_dig); 2141 for (; absValue >= 1; absValue *= T(0.5)) 2142 expCount++; 2143 } 2144 exp = expCount; 2145 return value < 0 ? -absValue : absValue; 2146 } 2147 2148 Unqual!T vf = value; 2149 ushort* vu = cast(ushort*)&vf; 2150 static if (is(immutable T == immutable float)) 2151 int* vi = cast(int*)&vf; 2152 else 2153 long* vl = cast(long*)&vf; 2154 int ex; 2155 alias F = floatTraits!T; 2156 2157 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2158 static if (F.realFormat == RealFormat.ieeeExtended || 2159 F.realFormat == RealFormat.ieeeExtended53) 2160 { 2161 if (ex) 2162 { // If exponent is non-zero 2163 if (ex == F.EXPMASK) // infinity or NaN 2164 { 2165 if (*vl & 0x7FFF_FFFF_FFFF_FFFF) // NaN 2166 { 2167 *vl |= 0xC000_0000_0000_0000; // convert NaNS to NaNQ 2168 exp = int.min; 2169 } 2170 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity 2171 exp = int.min; 2172 else // positive infinity 2173 exp = int.max; 2174 2175 } 2176 else 2177 { 2178 exp = ex - F.EXPBIAS; 2179 vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE; 2180 } 2181 } 2182 else if (!*vl) 2183 { 2184 // vf is +-0.0 2185 exp = 0; 2186 } 2187 else 2188 { 2189 // subnormal 2190 2191 vf *= F.RECIP_EPSILON; 2192 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2193 exp = ex - F.EXPBIAS - T.mant_dig + 1; 2194 vu[F.EXPPOS_SHORT] = ((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FFE; 2195 } 2196 return vf; 2197 } 2198 else static if (F.realFormat == RealFormat.ieeeQuadruple) 2199 { 2200 if (ex) // If exponent is non-zero 2201 { 2202 if (ex == F.EXPMASK) 2203 { 2204 // infinity or NaN 2205 if (vl[MANTISSA_LSB] | 2206 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN 2207 { 2208 // convert NaNS to NaNQ 2209 vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000; 2210 exp = int.min; 2211 } 2212 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity 2213 exp = int.min; 2214 else // positive infinity 2215 exp = int.max; 2216 } 2217 else 2218 { 2219 exp = ex - F.EXPBIAS; 2220 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]); 2221 } 2222 } 2223 else if ((vl[MANTISSA_LSB] | 2224 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) 2225 { 2226 // vf is +-0.0 2227 exp = 0; 2228 } 2229 else 2230 { 2231 // subnormal 2232 vf *= F.RECIP_EPSILON; 2233 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2234 exp = ex - F.EXPBIAS - T.mant_dig + 1; 2235 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]); 2236 } 2237 return vf; 2238 } 2239 else static if (F.realFormat == RealFormat.ieeeDouble) 2240 { 2241 if (ex) // If exponent is non-zero 2242 { 2243 if (ex == F.EXPMASK) // infinity or NaN 2244 { 2245 if (*vl == 0x7FF0_0000_0000_0000) // positive infinity 2246 { 2247 exp = int.max; 2248 } 2249 else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity 2250 exp = int.min; 2251 else 2252 { // NaN 2253 *vl |= 0x0008_0000_0000_0000; // convert NaNS to NaNQ 2254 exp = int.min; 2255 } 2256 } 2257 else 2258 { 2259 exp = (ex - F.EXPBIAS) >> 4; 2260 vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0); 2261 } 2262 } 2263 else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF)) 2264 { 2265 // vf is +-0.0 2266 exp = 0; 2267 } 2268 else 2269 { 2270 // subnormal 2271 vf *= F.RECIP_EPSILON; 2272 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2273 exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1; 2274 vu[F.EXPPOS_SHORT] = 2275 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FE0); 2276 } 2277 return vf; 2278 } 2279 else static if (F.realFormat == RealFormat.ieeeSingle) 2280 { 2281 if (ex) // If exponent is non-zero 2282 { 2283 if (ex == F.EXPMASK) // infinity or NaN 2284 { 2285 if (*vi == 0x7F80_0000) // positive infinity 2286 { 2287 exp = int.max; 2288 } 2289 else if (*vi == 0xFF80_0000) // negative infinity 2290 exp = int.min; 2291 else 2292 { // NaN 2293 *vi |= 0x0040_0000; // convert NaNS to NaNQ 2294 exp = int.min; 2295 } 2296 } 2297 else 2298 { 2299 exp = (ex - F.EXPBIAS) >> 7; 2300 vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00); 2301 } 2302 } 2303 else if (!(*vi & 0x7FFF_FFFF)) 2304 { 2305 // vf is +-0.0 2306 exp = 0; 2307 } 2308 else 2309 { 2310 // subnormal 2311 vf *= F.RECIP_EPSILON; 2312 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 2313 exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1; 2314 vu[F.EXPPOS_SHORT] = 2315 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00); 2316 } 2317 return vf; 2318 } 2319 else // static if (F.realFormat == RealFormat.ibmExtended) 2320 { 2321 assert(0, "frexp not implemented"); 2322 } 2323 } 2324 2325 /// 2326 @safe unittest 2327 { 2328 import std.math.operations : isClose; 2329 2330 int exp; 2331 real mantissa = frexp(123.456L, exp); 2332 2333 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L)); 2334 2335 assert(frexp(-real.nan, exp) && exp == int.min); 2336 assert(frexp(real.nan, exp) && exp == int.min); 2337 assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min); 2338 assert(frexp(real.infinity, exp) == real.infinity && exp == int.max); 2339 assert(frexp(-0.0, exp) == -0.0 && exp == 0); 2340 assert(frexp(0.0, exp) == 0.0 && exp == 0); 2341 } 2342 2343 @safe @nogc nothrow unittest 2344 { 2345 import std.math.operations : isClose; 2346 2347 int exp; 2348 real mantissa = frexp(123.456L, exp); 2349 2350 // check if values are equal to 19 decimal digits of precision 2351 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L, 1e-18)); 2352 } 2353 2354 @safe unittest 2355 { 2356 import std.math : floatTraits, RealFormat; 2357 import std.math.traits : isIdentical; 2358 import std.meta : AliasSeq; 2359 import std.typecons : tuple, Tuple; 2360 2361 static foreach (T; AliasSeq!(real, double, float)) 2362 {{ 2363 Tuple!(T, T, int)[] vals = [ // x,frexp,exp 2364 tuple(T(0.0), T( 0.0 ), 0), 2365 tuple(T(-0.0), T( -0.0), 0), 2366 tuple(T(1.0), T( .5 ), 1), 2367 tuple(T(-1.0), T( -.5 ), 1), 2368 tuple(T(2.0), T( .5 ), 2), 2369 tuple(T(float.min_normal/2.0f), T(.5), -126), 2370 tuple(T.infinity, T.infinity, int.max), 2371 tuple(-T.infinity, -T.infinity, int.min), 2372 tuple(T.nan, T.nan, int.min), 2373 tuple(-T.nan, -T.nan, int.min), 2374 2375 // https://issues.dlang.org/show_bug.cgi?id=16026: 2376 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2) 2377 ]; 2378 2379 foreach (elem; vals) 2380 { 2381 T x = elem[0]; 2382 T e = elem[1]; 2383 int exp = elem[2]; 2384 int eptr; 2385 T v = frexp(x, eptr); 2386 assert(isIdentical(e, v)); 2387 assert(exp == eptr); 2388 } 2389 2390 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended) 2391 { 2392 static T[3][] extendedvals = [ // x,frexp,exp 2393 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal 2394 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063], 2395 [T.min_normal, .5, -16381], 2396 [T.min_normal/2.0L, .5, -16382] // subnormal 2397 ]; 2398 foreach (elem; extendedvals) 2399 { 2400 T x = elem[0]; 2401 T e = elem[1]; 2402 int exp = cast(int) elem[2]; 2403 int eptr; 2404 T v = frexp(x, eptr); 2405 assert(isIdentical(e, v)); 2406 assert(exp == eptr); 2407 } 2408 } 2409 }} 2410 2411 // CTFE 2412 alias CtfeFrexpResult= Tuple!(real, int); 2413 static CtfeFrexpResult ctfeFrexp(T)(const T value) 2414 { 2415 int exp; 2416 auto significand = frexp(value, exp); 2417 return CtfeFrexpResult(significand, exp); 2418 } 2419 static foreach (T; AliasSeq!(real, double, float)) 2420 {{ 2421 enum Tuple!(T, T, int)[] vals = [ // x,frexp,exp 2422 tuple(T(0.0), T( 0.0 ), 0), 2423 tuple(T(-0.0), T( -0.0), 0), 2424 tuple(T(1.0), T( .5 ), 1), 2425 tuple(T(-1.0), T( -.5 ), 1), 2426 tuple(T(2.0), T( .5 ), 2), 2427 tuple(T(float.min_normal/2.0f), T(.5), -126), 2428 tuple(T.infinity, T.infinity, int.max), 2429 tuple(-T.infinity, -T.infinity, int.min), 2430 tuple(T.nan, T.nan, int.min), 2431 tuple(-T.nan, -T.nan, int.min), 2432 2433 // https://issues.dlang.org/show_bug.cgi?id=16026: 2434 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2) 2435 ]; 2436 2437 static foreach (elem; vals) 2438 { 2439 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], elem[2])); 2440 } 2441 2442 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended) 2443 { 2444 enum T[3][] extendedvals = [ // x,frexp,exp 2445 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal 2446 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063], 2447 [T.min_normal, .5, -16381], 2448 [T.min_normal/2.0L, .5, -16382] // subnormal 2449 ]; 2450 static foreach (elem; extendedvals) 2451 { 2452 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], cast(int) elem[2])); 2453 } 2454 } 2455 }} 2456 } 2457 2458 @safe unittest 2459 { 2460 import std.meta : AliasSeq; 2461 void foo() { 2462 static foreach (T; AliasSeq!(real, double, float)) 2463 {{ 2464 int exp; 2465 const T a = 1; 2466 immutable T b = 2; 2467 auto c = frexp(a, exp); 2468 auto d = frexp(b, exp); 2469 }} 2470 } 2471 } 2472 2473 /****************************************** 2474 * Extracts the exponent of x as a signed integral value. 2475 * 2476 * If x is not a special value, the result is the same as 2477 * $(D cast(int) logb(x)). 2478 * 2479 * $(TABLE_SV 2480 * $(TR $(TH x) $(TH ilogb(x)) $(TH Range error?)) 2481 * $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes)) 2482 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max) $(TD no)) 2483 * $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD no)) 2484 * ) 2485 */ 2486 int ilogb(T)(const T x) @trusted pure nothrow @nogc 2487 if (isFloatingPoint!T) 2488 { 2489 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB; 2490 2491 import core.bitop : bsr; 2492 alias F = floatTraits!T; 2493 2494 union floatBits 2495 { 2496 T rv; 2497 ushort[T.sizeof/2] vu; 2498 uint[T.sizeof/4] vui; 2499 static if (T.sizeof >= 8) 2500 ulong[T.sizeof/8] vul; 2501 } 2502 floatBits y = void; 2503 y.rv = x; 2504 2505 int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK; 2506 static if (F.realFormat == RealFormat.ieeeExtended || 2507 F.realFormat == RealFormat.ieeeExtended53) 2508 { 2509 if (ex) 2510 { 2511 // If exponent is non-zero 2512 if (ex == F.EXPMASK) // infinity or NaN 2513 { 2514 if (y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) // NaN 2515 return FP_ILOGBNAN; 2516 else // +-infinity 2517 return int.max; 2518 } 2519 else 2520 { 2521 return ex - F.EXPBIAS - 1; 2522 } 2523 } 2524 else if (!y.vul[0]) 2525 { 2526 // vf is +-0.0 2527 return FP_ILOGB0; 2528 } 2529 else 2530 { 2531 // subnormal 2532 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]); 2533 } 2534 } 2535 else static if (F.realFormat == RealFormat.ieeeQuadruple) 2536 { 2537 if (ex) // If exponent is non-zero 2538 { 2539 if (ex == F.EXPMASK) 2540 { 2541 // infinity or NaN 2542 if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN 2543 return FP_ILOGBNAN; 2544 else // +- infinity 2545 return int.max; 2546 } 2547 else 2548 { 2549 return ex - F.EXPBIAS - 1; 2550 } 2551 } 2552 else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) 2553 { 2554 // vf is +-0.0 2555 return FP_ILOGB0; 2556 } 2557 else 2558 { 2559 // subnormal 2560 const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF; 2561 const ulong lsb = y.vul[MANTISSA_LSB]; 2562 if (msb) 2563 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64; 2564 else 2565 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb); 2566 } 2567 } 2568 else static if (F.realFormat == RealFormat.ieeeDouble) 2569 { 2570 if (ex) // If exponent is non-zero 2571 { 2572 if (ex == F.EXPMASK) // infinity or NaN 2573 { 2574 if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000) // +- infinity 2575 return int.max; 2576 else // NaN 2577 return FP_ILOGBNAN; 2578 } 2579 else 2580 { 2581 return ((ex - F.EXPBIAS) >> 4) - 1; 2582 } 2583 } 2584 else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF)) 2585 { 2586 // vf is +-0.0 2587 return FP_ILOGB0; 2588 } 2589 else 2590 { 2591 // subnormal 2592 enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF; 2593 return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64); 2594 } 2595 } 2596 else static if (F.realFormat == RealFormat.ieeeSingle) 2597 { 2598 if (ex) // If exponent is non-zero 2599 { 2600 if (ex == F.EXPMASK) // infinity or NaN 2601 { 2602 if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000) // +- infinity 2603 return int.max; 2604 else // NaN 2605 return FP_ILOGBNAN; 2606 } 2607 else 2608 { 2609 return ((ex - F.EXPBIAS) >> 7) - 1; 2610 } 2611 } 2612 else if (!(y.vui[0] & 0x7FFF_FFFF)) 2613 { 2614 // vf is +-0.0 2615 return FP_ILOGB0; 2616 } 2617 else 2618 { 2619 // subnormal 2620 const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT; 2621 return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa); 2622 } 2623 } 2624 else // static if (F.realFormat == RealFormat.ibmExtended) 2625 { 2626 assert(0, "ilogb not implemented"); 2627 } 2628 } 2629 /// ditto 2630 int ilogb(T)(const T x) @safe pure nothrow @nogc 2631 if (isIntegral!T && isUnsigned!T) 2632 { 2633 import core.bitop : bsr; 2634 if (x == 0) 2635 return FP_ILOGB0; 2636 else 2637 { 2638 static assert(T.sizeof <= ulong.sizeof, "integer size too large for the current ilogb implementation"); 2639 return bsr(x); 2640 } 2641 } 2642 /// ditto 2643 int ilogb(T)(const T x) @safe pure nothrow @nogc 2644 if (isIntegral!T && isSigned!T) 2645 { 2646 import std.traits : Unsigned; 2647 // Note: abs(x) can not be used because the return type is not Unsigned and 2648 // the return value would be wrong for x == int.min 2649 Unsigned!T absx = x >= 0 ? x : -x; 2650 return ilogb(absx); 2651 } 2652 2653 /// 2654 @safe pure unittest 2655 { 2656 assert(ilogb(1) == 0); 2657 assert(ilogb(3) == 1); 2658 assert(ilogb(3.0) == 1); 2659 assert(ilogb(100_000_000) == 26); 2660 2661 assert(ilogb(0) == FP_ILOGB0); 2662 assert(ilogb(0.0) == FP_ILOGB0); 2663 assert(ilogb(double.nan) == FP_ILOGBNAN); 2664 assert(ilogb(double.infinity) == int.max); 2665 } 2666 2667 /** 2668 Special return values of $(LREF ilogb). 2669 */ 2670 alias FP_ILOGB0 = core.stdc.math.FP_ILOGB0; 2671 /// ditto 2672 alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN; 2673 2674 /// 2675 @safe pure unittest 2676 { 2677 assert(ilogb(0) == FP_ILOGB0); 2678 assert(ilogb(0.0) == FP_ILOGB0); 2679 assert(ilogb(double.nan) == FP_ILOGBNAN); 2680 } 2681 2682 @safe nothrow @nogc unittest 2683 { 2684 import std.math : floatTraits, RealFormat; 2685 import std.math.operations : nextUp; 2686 import std.meta : AliasSeq; 2687 import std.typecons : Tuple; 2688 static foreach (F; AliasSeq!(float, double, real)) 2689 {{ 2690 alias T = Tuple!(F, int); 2691 T[13] vals = // x, ilogb(x) 2692 [ 2693 T( F.nan , FP_ILOGBNAN ), 2694 T( -F.nan , FP_ILOGBNAN ), 2695 T( F.infinity, int.max ), 2696 T( -F.infinity, int.max ), 2697 T( 0.0 , FP_ILOGB0 ), 2698 T( -0.0 , FP_ILOGB0 ), 2699 T( 2.0 , 1 ), 2700 T( 2.0001 , 1 ), 2701 T( 1.9999 , 0 ), 2702 T( 0.5 , -1 ), 2703 T( 123.123 , 6 ), 2704 T( -123.123 , 6 ), 2705 T( 0.123 , -4 ), 2706 ]; 2707 2708 foreach (elem; vals) 2709 { 2710 assert(ilogb(elem[0]) == elem[1]); 2711 } 2712 }} 2713 2714 // min_normal and subnormals 2715 assert(ilogb(-float.min_normal) == -126); 2716 assert(ilogb(nextUp(-float.min_normal)) == -127); 2717 assert(ilogb(nextUp(-float(0.0))) == -149); 2718 assert(ilogb(-double.min_normal) == -1022); 2719 assert(ilogb(nextUp(-double.min_normal)) == -1023); 2720 assert(ilogb(nextUp(-double(0.0))) == -1074); 2721 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended) 2722 { 2723 assert(ilogb(-real.min_normal) == -16382); 2724 assert(ilogb(nextUp(-real.min_normal)) == -16383); 2725 assert(ilogb(nextUp(-real(0.0))) == -16445); 2726 } 2727 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) 2728 { 2729 assert(ilogb(-real.min_normal) == -1022); 2730 assert(ilogb(nextUp(-real.min_normal)) == -1023); 2731 assert(ilogb(nextUp(-real(0.0))) == -1074); 2732 } 2733 2734 // test integer types 2735 assert(ilogb(0) == FP_ILOGB0); 2736 assert(ilogb(int.max) == 30); 2737 assert(ilogb(int.min) == 31); 2738 assert(ilogb(uint.max) == 31); 2739 assert(ilogb(long.max) == 62); 2740 assert(ilogb(long.min) == 63); 2741 assert(ilogb(ulong.max) == 63); 2742 } 2743 2744 /******************************************* 2745 * Compute n * 2$(SUPERSCRIPT exp) 2746 * References: frexp 2747 */ 2748 2749 pragma(inline, true) 2750 real ldexp(real n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); } 2751 ///ditto 2752 pragma(inline, true) 2753 double ldexp(double n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); } 2754 ///ditto 2755 pragma(inline, true) 2756 float ldexp(float n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); } 2757 2758 /// 2759 @nogc @safe pure nothrow unittest 2760 { 2761 import std.meta : AliasSeq; 2762 static foreach (T; AliasSeq!(float, double, real)) 2763 {{ 2764 T r; 2765 2766 r = ldexp(3.0L, 3); 2767 assert(r == 24); 2768 2769 r = ldexp(cast(T) 3.0, cast(int) 3); 2770 assert(r == 24); 2771 2772 T n = 3.0; 2773 int exp = 3; 2774 r = ldexp(n, exp); 2775 assert(r == 24); 2776 }} 2777 } 2778 2779 @safe pure nothrow @nogc unittest 2780 { 2781 import std.math : floatTraits, RealFormat; 2782 2783 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended || 2784 floatTraits!(real).realFormat == RealFormat.ieeeExtended53 || 2785 floatTraits!(real).realFormat == RealFormat.ieeeQuadruple) 2786 { 2787 assert(ldexp(1.0L, -16384) == 0x1p-16384L); 2788 assert(ldexp(1.0L, -16382) == 0x1p-16382L); 2789 int x; 2790 real n = frexp(0x1p-16384L, x); 2791 assert(n == 0.5L); 2792 assert(x==-16383); 2793 assert(ldexp(n, x)==0x1p-16384L); 2794 } 2795 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) 2796 { 2797 assert(ldexp(1.0L, -1024) == 0x1p-1024L); 2798 assert(ldexp(1.0L, -1022) == 0x1p-1022L); 2799 int x; 2800 real n = frexp(0x1p-1024L, x); 2801 assert(n == 0.5L); 2802 assert(x==-1023); 2803 assert(ldexp(n, x)==0x1p-1024L); 2804 } 2805 else static assert(false, "Floating point type real not supported"); 2806 } 2807 2808 /* workaround https://issues.dlang.org/show_bug.cgi?id=14718 2809 float parsing depends on platform strtold 2810 @safe pure nothrow @nogc unittest 2811 { 2812 assert(ldexp(1.0, -1024) == 0x1p-1024); 2813 assert(ldexp(1.0, -1022) == 0x1p-1022); 2814 int x; 2815 double n = frexp(0x1p-1024, x); 2816 assert(n == 0.5); 2817 assert(x==-1023); 2818 assert(ldexp(n, x)==0x1p-1024); 2819 } 2820 2821 @safe pure nothrow @nogc unittest 2822 { 2823 assert(ldexp(1.0f, -128) == 0x1p-128f); 2824 assert(ldexp(1.0f, -126) == 0x1p-126f); 2825 int x; 2826 float n = frexp(0x1p-128f, x); 2827 assert(n == 0.5f); 2828 assert(x==-127); 2829 assert(ldexp(n, x)==0x1p-128f); 2830 } 2831 */ 2832 2833 @safe @nogc nothrow unittest 2834 { 2835 import std.math.operations : isClose; 2836 2837 static real[3][] vals = // value,exp,ldexp 2838 [ 2839 [ 0, 0, 0], 2840 [ 1, 0, 1], 2841 [ -1, 0, -1], 2842 [ 1, 1, 2], 2843 [ 123, 10, 125952], 2844 [ real.max, int.max, real.infinity], 2845 [ real.max, -int.max, 0], 2846 [ real.min_normal, -int.max, 0], 2847 ]; 2848 int i; 2849 2850 for (i = 0; i < vals.length; i++) 2851 { 2852 real x = vals[i][0]; 2853 int exp = cast(int) vals[i][1]; 2854 real z = vals[i][2]; 2855 real l = ldexp(x, exp); 2856 2857 assert(isClose(z, l, 1e-6)); 2858 } 2859 2860 real function(real, int) pldexp = &ldexp; 2861 assert(pldexp != null); 2862 } 2863 2864 private 2865 { 2866 // Coefficients shared across log(), log2(), log10(), log1p(). 2867 template LogCoeffs(T) 2868 { 2869 import std.math : floatTraits, RealFormat; 2870 2871 static if (floatTraits!T.realFormat == RealFormat.ieeeQuadruple) 2872 { 2873 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x) 2874 // Theoretical peak relative error = 5.3e-37 2875 static immutable real[13] logP = [ 2876 1.313572404063446165910279910527789794488E4L, 2877 7.771154681358524243729929227226708890930E4L, 2878 2.014652742082537582487669938141683759923E5L, 2879 3.007007295140399532324943111654767187848E5L, 2880 2.854829159639697837788887080758954924001E5L, 2881 1.797628303815655343403735250238293741397E5L, 2882 7.594356839258970405033155585486712125861E4L, 2883 2.128857716871515081352991964243375186031E4L, 2884 3.824952356185897735160588078446136783779E3L, 2885 4.114517881637811823002128927449878962058E2L, 2886 2.321125933898420063925789532045674660756E1L, 2887 4.998469661968096229986658302195402690910E-1L, 2888 1.538612243596254322971797716843006400388E-6L 2889 ]; 2890 static immutable real[13] logQ = [ 2891 3.940717212190338497730839731583397586124E4L, 2892 2.626900195321832660448791748036714883242E5L, 2893 7.777690340007566932935753241556479363645E5L, 2894 1.347518538384329112529391120390701166528E6L, 2895 1.514882452993549494932585972882995548426E6L, 2896 1.158019977462989115839826904108208787040E6L, 2897 6.132189329546557743179177159925690841200E5L, 2898 2.248234257620569139969141618556349415120E5L, 2899 5.605842085972455027590989944010492125825E4L, 2900 9.147150349299596453976674231612674085381E3L, 2901 9.104928120962988414618126155557301584078E2L, 2902 4.839208193348159620282142911143429644326E1L, 2903 1.0 2904 ]; 2905 2906 // log2 uses the same coefficients as log. 2907 alias log2P = logP; 2908 alias log2Q = logQ; 2909 2910 // log10 uses the same coefficients as log. 2911 alias log10P = logP; 2912 alias log10Q = logQ; 2913 2914 // Coefficients for log(x) = z + z^^3 R(z^^2)/S(z^^2) 2915 // where z = 2(x-1)/(x+1) 2916 // Theoretical peak relative error = 1.1e-35 2917 static immutable real[6] logR = [ 2918 1.418134209872192732479751274970992665513E5L, 2919 -8.977257995689735303686582344659576526998E4L, 2920 2.048819892795278657810231591630928516206E4L, 2921 -2.024301798136027039250415126250455056397E3L, 2922 8.057002716646055371965756206836056074715E1L, 2923 -8.828896441624934385266096344596648080902E-1L 2924 ]; 2925 static immutable real[7] logS = [ 2926 1.701761051846631278975701529965589676574E6L, 2927 -1.332535117259762928288745111081235577029E6L, 2928 4.001557694070773974936904547424676279307E5L, 2929 -5.748542087379434595104154610899551484314E4L, 2930 3.998526750980007367835804959888064681098E3L, 2931 -1.186359407982897997337150403816839480438E2L, 2932 1.0 2933 ]; 2934 } 2935 else static if (floatTraits!T.realFormat == RealFormat.ieeeExtended || 2936 floatTraits!T.realFormat == RealFormat.ieeeExtended53) 2937 { 2938 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x) 2939 // Theoretical peak relative error = 2.32e-20 2940 static immutable real[7] logP = [ 2941 2.0039553499201281259648E1L, 2942 5.7112963590585538103336E1L, 2943 6.0949667980987787057556E1L, 2944 2.9911919328553073277375E1L, 2945 6.5787325942061044846969E0L, 2946 4.9854102823193375972212E-1L, 2947 4.5270000862445199635215E-5L, 2948 ]; 2949 static immutable real[7] logQ = [ 2950 6.0118660497603843919306E1L, 2951 2.1642788614495947685003E2L, 2952 3.0909872225312059774938E2L, 2953 2.2176239823732856465394E2L, 2954 8.3047565967967209469434E1L, 2955 1.5062909083469192043167E1L, 2956 1.0000000000000000000000E0L, 2957 ]; 2958 2959 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x) 2960 // Theoretical peak relative error = 6.2e-22 2961 static immutable real[7] log2P = [ 2962 1.0747524399916215149070E2L, 2963 3.4258224542413922935104E2L, 2964 4.2401812743503691187826E2L, 2965 2.5620629828144409632571E2L, 2966 7.7671073698359539859595E1L, 2967 1.0767376367209449010438E1L, 2968 4.9962495940332550844739E-1L, 2969 ]; 2970 static immutable real[8] log2Q = [ 2971 3.2242573199748645407652E2L, 2972 1.2695660352705325274404E3L, 2973 2.0307734695595183428202E3L, 2974 1.6911722418503949084863E3L, 2975 7.7952888181207260646090E2L, 2976 1.9444210022760132894510E2L, 2977 2.3479774160285863271658E1L, 2978 1.0000000000000000000000E0, 2979 ]; 2980 2981 // log10 uses the same coefficients as log2. 2982 alias log10P = log2P; 2983 alias log10Q = log2Q; 2984 2985 // Coefficients for log(x) = z + z^^3 R(z^^2)/S(z^^2) 2986 // where z = 2(x-1)/(x+1) 2987 // Theoretical peak relative error = 6.16e-22 2988 static immutable real[4] logR = [ 2989 -3.5717684488096787370998E1L, 2990 1.0777257190312272158094E1L, 2991 -7.1990767473014147232598E-1L, 2992 1.9757429581415468984296E-3L, 2993 ]; 2994 static immutable real[4] logS = [ 2995 -4.2861221385716144629696E2L, 2996 1.9361891836232102174846E2L, 2997 -2.6201045551331104417768E1L, 2998 1.0000000000000000000000E0L, 2999 ]; 3000 } 3001 else static if (floatTraits!T.realFormat == RealFormat.ieeeDouble) 3002 { 3003 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x) 3004 static immutable double[6] logP = [ 3005 7.70838733755885391666E0, 3006 1.79368678507819816313E1, 3007 1.44989225341610930846E1, 3008 4.70579119878881725854E0, 3009 4.97494994976747001425E-1, 3010 1.01875663804580931796E-4, 3011 ]; 3012 static immutable double[6] logQ = [ 3013 2.31251620126765340583E1, 3014 7.11544750618563894466E1, 3015 8.29875266912776603211E1, 3016 4.52279145837532221105E1, 3017 1.12873587189167450590E1, 3018 1.00000000000000000000E0, 3019 ]; 3020 3021 // log2 uses the same coefficients as log. 3022 alias log2P = logP; 3023 alias log2Q = logQ; 3024 3025 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x) 3026 static immutable double[7] logp1P = [ 3027 2.0039553499201281259648E1, 3028 5.7112963590585538103336E1, 3029 6.0949667980987787057556E1, 3030 2.9911919328553073277375E1, 3031 6.5787325942061044846969E0, 3032 4.9854102823193375972212E-1, 3033 4.5270000862445199635215E-5, 3034 ]; 3035 static immutable double[7] logp1Q = [ 3036 6.0118660497603843919306E1, 3037 2.1642788614495947685003E2, 3038 3.0909872225312059774938E2, 3039 2.2176239823732856465394E2, 3040 8.3047565967967209469434E1, 3041 1.5062909083469192043167E1, 3042 1.0000000000000000000000E0, 3043 ]; 3044 3045 static immutable double[7] log10P = [ 3046 1.98892446572874072159E1, 3047 5.67349287391754285487E1, 3048 6.06127134467767258030E1, 3049 2.97877425097986925891E1, 3050 6.56312093769992875930E0, 3051 4.98531067254050724270E-1, 3052 4.58482948458143443514E-5, 3053 ]; 3054 static immutable double[7] log10Q = [ 3055 5.96677339718622216300E1, 3056 2.14955586696422947765E2, 3057 3.07254189979530058263E2, 3058 2.20664384982121929218E2, 3059 8.27410449222435217021E1, 3060 1.50314182634250003249E1, 3061 1.00000000000000000000E0, 3062 ]; 3063 3064 // Coefficients for log(x) = z + z^^3 R(z)/S(z) 3065 // where z = 2(x-1)/(x+1) 3066 static immutable double[3] logR = [ 3067 -6.41409952958715622951E1, 3068 1.63866645699558079767E1, 3069 -7.89580278884799154124E-1, 3070 ]; 3071 static immutable double[4] logS = [ 3072 -7.69691943550460008604E2, 3073 3.12093766372244180303E2, 3074 -3.56722798256324312549E1, 3075 1.00000000000000000000E0, 3076 ]; 3077 } 3078 else static if (floatTraits!T.realFormat == RealFormat.ieeeSingle) 3079 { 3080 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x) 3081 static immutable float[9] logP = [ 3082 3.3333331174E-1, 3083 -2.4999993993E-1, 3084 2.0000714765E-1, 3085 -1.6668057665E-1, 3086 1.4249322787E-1, 3087 -1.2420140846E-1, 3088 1.1676998740E-1, 3089 -1.1514610310E-1, 3090 7.0376836292E-2, 3091 ]; 3092 3093 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x) 3094 static immutable float[7] logp1P = [ 3095 2.0039553499E1, 3096 5.7112963590E1, 3097 6.0949667980E1, 3098 2.9911919328E1, 3099 6.5787325942E0, 3100 4.9854102823E-1, 3101 4.5270000862E-5, 3102 ]; 3103 static immutable float[7] logp1Q = [ 3104 6.01186604976E1, 3105 2.16427886144E2, 3106 3.09098722253E2, 3107 2.21762398237E2, 3108 8.30475659679E1, 3109 1.50629090834E1, 3110 1.00000000000E0, 3111 ]; 3112 3113 // log2 and log10 uses the same coefficients as log. 3114 alias log2P = logP; 3115 alias log10P = logP; 3116 } 3117 else 3118 static assert(0, "no coefficients for log()"); 3119 } 3120 } 3121 3122 /************************************** 3123 * Calculate the natural logarithm of x. 3124 * 3125 * $(TABLE_SV 3126 * $(TR $(TH x) $(TH log(x)) $(TH divide by 0?) $(TH invalid?)) 3127 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) 3128 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes)) 3129 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) 3130 * ) 3131 */ 3132 pragma(inline, true) 3133 real log(real x) @safe pure nothrow @nogc 3134 { 3135 version (INLINE_YL2X) 3136 { 3137 import std.math.constants : LN2; 3138 return core.math.yl2x(x, LN2); 3139 } 3140 else 3141 return logImpl(x); 3142 } 3143 3144 /// ditto 3145 pragma(inline, true) 3146 double log(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log(cast(real) x) : logImpl(x); } 3147 3148 /// ditto 3149 pragma(inline, true) 3150 float log(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log(cast(real) x) : logImpl(x); } 3151 3152 // @@@DEPRECATED_[2.112.0]@@@ 3153 deprecated("`std.math.exponential.log` called with argument types `(int)` matches both " 3154 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.") 3155 real log(int x) @safe pure nothrow @nogc { return log(cast(real) x); } 3156 // @@@DEPRECATED_[2.112.0]@@@ 3157 deprecated("`std.math.exponential.log` called with argument types `(uint)` matches both " 3158 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.") 3159 real log(uint x) @safe pure nothrow @nogc { return log(cast(real) x); } 3160 // @@@DEPRECATED_[2.112.0]@@@ 3161 deprecated("`std.math.exponential.log` called with argument types `(long)` matches both " 3162 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.") 3163 real log(long x) @safe pure nothrow @nogc { return log(cast(real) x); } 3164 // @@@DEPRECATED_[2.112.0]@@@ 3165 deprecated("`std.math.exponential.log` called with argument types `(ulong)` matches both " 3166 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.") 3167 real log(ulong x) @safe pure nothrow @nogc { return log(cast(real) x); } 3168 3169 /// 3170 @safe pure nothrow @nogc unittest 3171 { 3172 import std.math.operations : feqrel; 3173 import std.math.constants : E; 3174 3175 assert(feqrel(log(E), 1) >= real.mant_dig - 1); 3176 } 3177 3178 private T logImpl(T, bool LOG1P = false)(T x) @safe pure nothrow @nogc 3179 { 3180 import std.math.constants : SQRT1_2; 3181 import std.math.algebraic : poly; 3182 import std.math.traits : isInfinity, isNaN, signbit; 3183 import std.math : floatTraits, RealFormat; 3184 3185 alias coeffs = LogCoeffs!T; 3186 alias F = floatTraits!T; 3187 3188 static if (LOG1P) 3189 { 3190 const T xm1 = x; 3191 x = x + 1.0; 3192 } 3193 3194 static if (F.realFormat == RealFormat.ieeeExtended || 3195 F.realFormat == RealFormat.ieeeExtended53 || 3196 F.realFormat == RealFormat.ieeeQuadruple) 3197 { 3198 // C1 + C2 = LN2. 3199 enum T C1 = 6.93145751953125E-1L; 3200 enum T C2 = 1.428606820309417232121458176568075500134E-6L; 3201 } 3202 else static if (F.realFormat == RealFormat.ieeeDouble) 3203 { 3204 enum T C1 = 0.693359375; 3205 enum T C2 = -2.121944400546905827679e-4; 3206 } 3207 else static if (F.realFormat == RealFormat.ieeeSingle) 3208 { 3209 enum T C1 = 0.693359375; 3210 enum T C2 = -2.12194440e-4; 3211 } 3212 else 3213 static assert(0, "Not implemented for this architecture"); 3214 3215 // Special cases. 3216 if (isNaN(x)) 3217 return x; 3218 if (isInfinity(x) && !signbit(x)) 3219 return x; 3220 if (x == 0.0) 3221 return -T.infinity; 3222 if (x < 0.0) 3223 return T.nan; 3224 3225 // Separate mantissa from exponent. 3226 // Note, frexp is used so that denormal numbers will be handled properly. 3227 T y, z; 3228 int exp; 3229 3230 x = frexp(x, exp); 3231 3232 static if (F.realFormat == RealFormat.ieeeDouble || 3233 F.realFormat == RealFormat.ieeeExtended || 3234 F.realFormat == RealFormat.ieeeExtended53 || 3235 F.realFormat == RealFormat.ieeeQuadruple) 3236 { 3237 // Logarithm using log(x) = z + z^^3 R(z) / S(z), 3238 // where z = 2(x - 1)/(x + 1) 3239 if ((exp > 2) || (exp < -2)) 3240 { 3241 if (x < SQRT1_2) 3242 { // 2(2x - 1)/(2x + 1) 3243 exp -= 1; 3244 z = x - 0.5; 3245 y = 0.5 * z + 0.5; 3246 } 3247 else 3248 { // 2(x - 1)/(x + 1) 3249 z = x - 0.5; 3250 z -= 0.5; 3251 y = 0.5 * x + 0.5; 3252 } 3253 x = z / y; 3254 z = x * x; 3255 z = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS)); 3256 z += exp * C2; 3257 z += x; 3258 z += exp * C1; 3259 3260 return z; 3261 } 3262 } 3263 3264 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x) 3265 if (x < SQRT1_2) 3266 { 3267 exp -= 1; 3268 static if (LOG1P) 3269 { 3270 if (exp != 0) 3271 x = 2.0 * x - 1.0; 3272 else 3273 x = xm1; 3274 } 3275 else 3276 x = 2.0 * x - 1.0; 3277 3278 } 3279 else 3280 { 3281 static if (LOG1P) 3282 { 3283 if (exp != 0) 3284 x = x - 1.0; 3285 else 3286 x = xm1; 3287 } 3288 else 3289 x = x - 1.0; 3290 } 3291 z = x * x; 3292 static if (F.realFormat == RealFormat.ieeeSingle) 3293 y = x * (z * poly(x, coeffs.logP)); 3294 else 3295 y = x * (z * poly(x, coeffs.logP) / poly(x, coeffs.logQ)); 3296 y += exp * C2; 3297 z = y - 0.5 * z; 3298 3299 // Note, the sum of above terms does not exceed x/4, 3300 // so it contributes at most about 1/4 lsb to the error. 3301 z += x; 3302 z += exp * C1; 3303 3304 return z; 3305 } 3306 3307 @safe @nogc nothrow unittest 3308 { 3309 import std.math : floatTraits, RealFormat; 3310 import std.meta : AliasSeq; 3311 3312 static void testLog(T)(T[2][] vals) 3313 { 3314 import std.math.operations : isClose; 3315 import std.math.traits : isNaN; 3316 foreach (ref pair; vals) 3317 { 3318 if (isNaN(pair[1])) 3319 assert(isNaN(log(pair[0]))); 3320 else 3321 assert(isClose(log(pair[0]), pair[1])); 3322 } 3323 } 3324 static foreach (F; AliasSeq!(float, double, real)) 3325 {{ 3326 scope F[2][] vals = [ 3327 [F(1), F(0x0p+0)], [F(2), F(0x1.62e42fefa39ef358p-1)], 3328 [F(4), F(0x1.62e42fefa39ef358p+0)], [F(8), F(0x1.0a2b23f3bab73682p+1)], 3329 [F(16), F(0x1.62e42fefa39ef358p+1)], [F(32), F(0x1.bb9d3beb8c86b02ep+1)], 3330 [F(64), F(0x1.0a2b23f3bab73682p+2)], [F(128), F(0x1.3687a9f1af2b14ecp+2)], 3331 [F(256), F(0x1.62e42fefa39ef358p+2)], [F(512), F(0x1.8f40b5ed9812d1c2p+2)], 3332 [F(1024), F(0x1.bb9d3beb8c86b02ep+2)], [F(2048), F(0x1.e7f9c1e980fa8e98p+2)], 3333 [F(3), F(0x1.193ea7aad030a976p+0)], [F(5), F(0x1.9c041f7ed8d336bp+0)], 3334 [F(7), F(0x1.f2272ae325a57546p+0)], [F(15), F(0x1.5aa16394d481f014p+1)], 3335 [F(17), F(0x1.6aa6bc1fa7f79cfp+1)], [F(31), F(0x1.b78ce48912b59f12p+1)], 3336 [F(33), F(0x1.bf8d8f4d5b8d1038p+1)], [F(63), F(0x1.09291e8e3181b20ep+2)], 3337 [F(65), F(0x1.0b292939429755ap+2)], [F(-0), -F.infinity], [F(0), -F.infinity], 3338 [F(0.1), F(-0x1.26bb1bbb5551582ep+1)], [F(0.25), F(-0x1.62e42fefa39ef358p+0)], 3339 [F(0.75), F(-0x1.269621134db92784p-2)], [F(0.875), F(-0x1.1178e8227e47bde4p-3)], 3340 [F(10), F(0x1.26bb1bbb5551582ep+1)], [F(100), F(0x1.26bb1bbb5551582ep+2)], 3341 [F(10000), F(0x1.26bb1bbb5551582ep+3)], 3342 ]; 3343 testLog(vals); 3344 }} 3345 { 3346 float[2][16] vals = [ 3347 [float.nan, float.nan],[-float.nan, float.nan], 3348 [float.infinity, float.infinity], [-float.infinity, float.nan], 3349 [float.min_normal, -0x1.5d58ap+6f], [-float.min_normal, float.nan], 3350 [float.max, 0x1.62e43p+6f], [-float.max, float.nan], 3351 [float.min_normal / 2, -0x1.601e68p+6f], [-float.min_normal / 2, float.nan], 3352 [float.max / 2, 0x1.601e68p+6f], [-float.max / 2, float.nan], 3353 [float.min_normal / 3, -0x1.61bd9ap+6f], [-float.min_normal / 3, float.nan], 3354 [float.max / 3, 0x1.5e7f36p+6f], [-float.max / 3, float.nan], 3355 ]; 3356 testLog(vals); 3357 } 3358 { 3359 double[2][16] vals = [ 3360 [double.nan, double.nan],[-double.nan, double.nan], 3361 [double.infinity, double.infinity], [-double.infinity, double.nan], 3362 [double.min_normal, -0x1.6232bdd7abcd2p+9], [-double.min_normal, double.nan], 3363 [double.max, 0x1.62e42fefa39efp+9], [-double.max, double.nan], 3364 [double.min_normal / 2, -0x1.628b76e3a7b61p+9], [-double.min_normal / 2, double.nan], 3365 [double.max / 2, 0x1.628b76e3a7b61p+9], [-double.max / 2, double.nan], 3366 [double.min_normal / 3, -0x1.62bf5d2b81354p+9], [-double.min_normal / 3, double.nan], 3367 [double.max / 3, 0x1.6257909bce36ep+9], [-double.max / 3, double.nan], 3368 ]; 3369 testLog(vals); 3370 } 3371 alias F = floatTraits!real; 3372 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple) 3373 {{ 3374 real[2][16] vals = [ 3375 [real.nan, real.nan],[-real.nan, real.nan], 3376 [real.infinity, real.infinity], [-real.infinity, real.nan], 3377 [real.min_normal, -0x1.62d918ce2421d66p+13L], [-real.min_normal, real.nan], 3378 [real.max, 0x1.62e42fefa39ef358p+13L], [-real.max, real.nan], 3379 [real.min_normal / 2, -0x1.62dea45ee3e064dcp+13L], [-real.min_normal / 2, real.nan], 3380 [real.max / 2, 0x1.62dea45ee3e064dcp+13L], [-real.max / 2, real.nan], 3381 [real.min_normal / 3, -0x1.62e1e2c3617857e6p+13L], [-real.min_normal / 3, real.nan], 3382 [real.max / 3, 0x1.62db65fa664871d2p+13L], [-real.max / 3, real.nan], 3383 ]; 3384 testLog(vals); 3385 }} 3386 } 3387 3388 /************************************** 3389 * Calculate the base-10 logarithm of x. 3390 * 3391 * $(TABLE_SV 3392 * $(TR $(TH x) $(TH log10(x)) $(TH divide by 0?) $(TH invalid?)) 3393 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) 3394 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes)) 3395 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) 3396 * ) 3397 */ 3398 pragma(inline, true) 3399 real log10(real x) @safe pure nothrow @nogc 3400 { 3401 version (INLINE_YL2X) 3402 { 3403 import std.math.constants : LOG2; 3404 return core.math.yl2x(x, LOG2); 3405 } 3406 else 3407 return log10Impl(x); 3408 } 3409 3410 /// ditto 3411 pragma(inline, true) 3412 double log10(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log10(cast(real) x) : log10Impl(x); } 3413 3414 /// ditto 3415 pragma(inline, true) 3416 float log10(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log10(cast(real) x) : log10Impl(x); } 3417 3418 // @@@DEPRECATED_[2.112.0]@@@ 3419 deprecated("`std.math.exponential.log10` called with argument types `(int)` matches both " 3420 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.") 3421 real log10(int x) @safe pure nothrow @nogc { return log10(cast(real) x); } 3422 // @@@DEPRECATED_[2.112.0]@@@ 3423 deprecated("`std.math.exponential.log10` called with argument types `(uint)` matches both " 3424 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.") 3425 real log10(uint x) @safe pure nothrow @nogc { return log10(cast(real) x); } 3426 // @@@DEPRECATED_[2.112.0]@@@ 3427 deprecated("`std.math.exponential.log10` called with argument types `(long)` matches both " 3428 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.") 3429 real log10(long x) @safe pure nothrow @nogc { return log10(cast(real) x); } 3430 // @@@DEPRECATED_[2.112.0]@@@ 3431 deprecated("`std.math.exponential.log10` called with argument types `(ulong)` matches both " 3432 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.") 3433 real log10(ulong x) @safe pure nothrow @nogc { return log10(cast(real) x); } 3434 3435 /// 3436 @safe pure nothrow @nogc unittest 3437 { 3438 import std.math.algebraic : fabs; 3439 3440 assert(fabs(log10(1000.0L) - 3) < .000001); 3441 } 3442 3443 @safe pure nothrow @nogc unittest 3444 { 3445 import std.math.algebraic : fabs; 3446 3447 assert(fabs(log10(1000.0) - 3) < .000001); 3448 assert(fabs(log10(1000.0f) - 3) < .000001); 3449 } 3450 3451 private T log10Impl(T)(T x) @safe pure nothrow @nogc 3452 { 3453 import std.math.constants : SQRT1_2; 3454 import std.math.algebraic : poly; 3455 import std.math.traits : isNaN, isInfinity, signbit; 3456 import std.math : floatTraits, RealFormat; 3457 3458 alias coeffs = LogCoeffs!T; 3459 alias F = floatTraits!T; 3460 3461 static if (F.realFormat == RealFormat.ieeeExtended || 3462 F.realFormat == RealFormat.ieeeExtended53 || 3463 F.realFormat == RealFormat.ieeeQuadruple) 3464 { 3465 // log10(2) split into two parts. 3466 enum T L102A = 0.3125L; 3467 enum T L102B = -1.14700043360188047862611052755069732318101185E-2L; 3468 3469 // log10(e) split into two parts. 3470 enum T L10EA = 0.5L; 3471 enum T L10EB = -6.570551809674817234887108108339491770560299E-2L; 3472 } 3473 else static if (F.realFormat == RealFormat.ieeeDouble || 3474 F.realFormat == RealFormat.ieeeSingle) 3475 { 3476 enum T L102A = 3.0078125E-1; 3477 enum T L102B = 2.48745663981195213739E-4; 3478 3479 enum T L10EA = 4.3359375E-1; 3480 enum T L10EB = 7.00731903251827651129E-4; 3481 } 3482 else 3483 static assert(0, "Not implemented for this architecture"); 3484 3485 // Special cases are the same as for log. 3486 if (isNaN(x)) 3487 return x; 3488 if (isInfinity(x) && !signbit(x)) 3489 return x; 3490 if (x == 0.0) 3491 return -T.infinity; 3492 if (x < 0.0) 3493 return T.nan; 3494 3495 // Separate mantissa from exponent. 3496 // Note, frexp is used so that denormal numbers will be handled properly. 3497 T y, z; 3498 int exp; 3499 3500 x = frexp(x, exp); 3501 3502 static if (F.realFormat == RealFormat.ieeeExtended || 3503 F.realFormat == RealFormat.ieeeExtended53 || 3504 F.realFormat == RealFormat.ieeeQuadruple) 3505 { 3506 // Logarithm using log(x) = z + z^^3 R(z) / S(z), 3507 // where z = 2(x - 1)/(x + 1) 3508 if ((exp > 2) || (exp < -2)) 3509 { 3510 if (x < SQRT1_2) 3511 { // 2(2x - 1)/(2x + 1) 3512 exp -= 1; 3513 z = x - 0.5; 3514 y = 0.5 * z + 0.5; 3515 } 3516 else 3517 { // 2(x - 1)/(x + 1) 3518 z = x - 0.5; 3519 z -= 0.5; 3520 y = 0.5 * x + 0.5; 3521 } 3522 x = z / y; 3523 z = x * x; 3524 y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS)); 3525 goto Ldone; 3526 } 3527 } 3528 3529 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x) 3530 if (x < SQRT1_2) 3531 { 3532 exp -= 1; 3533 x = 2.0 * x - 1.0; 3534 } 3535 else 3536 x = x - 1.0; 3537 3538 z = x * x; 3539 static if (F.realFormat == RealFormat.ieeeSingle) 3540 y = x * (z * poly(x, coeffs.log10P)); 3541 else 3542 y = x * (z * poly(x, coeffs.log10P) / poly(x, coeffs.log10Q)); 3543 y = y - 0.5 * z; 3544 3545 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2). 3546 // This sequence of operations is critical and it may be horribly 3547 // defeated by some compiler optimizers. 3548 Ldone: 3549 z = y * L10EB; 3550 z += x * L10EB; 3551 z += exp * L102B; 3552 z += y * L10EA; 3553 z += x * L10EA; 3554 z += exp * L102A; 3555 3556 return z; 3557 } 3558 3559 @safe @nogc nothrow unittest 3560 { 3561 import std.math : floatTraits, RealFormat; 3562 import std.meta : AliasSeq; 3563 3564 static void testLog10(T)(T[2][] vals) 3565 { 3566 import std.math.operations : isClose; 3567 import std.math.traits : isNaN; 3568 foreach (ref pair; vals) 3569 { 3570 if (isNaN(pair[1])) 3571 assert(isNaN(log10(pair[0]))); 3572 else 3573 assert(isClose(log10(pair[0]), pair[1])); 3574 } 3575 } 3576 static foreach (F; AliasSeq!(float, double, real)) 3577 {{ 3578 scope F[2][] vals = [ 3579 [F(1), F(0x0p+0)], [F(2), F(0x1.34413509f79fef32p-2)], 3580 [F(4), F(0x1.34413509f79fef32p-1)], [F(8), F(0x1.ce61cf8ef36fe6cap-1)], 3581 [F(16), F(0x1.34413509f79fef32p+0)], [F(32), F(0x1.8151824c7587eafep+0)], 3582 [F(64), F(0x1.ce61cf8ef36fe6cap+0)], [F(128), F(0x1.0db90e68b8abf14cp+1)], 3583 [F(256), F(0x1.34413509f79fef32p+1)], [F(512), F(0x1.5ac95bab3693ed18p+1)], 3584 [F(1024), F(0x1.8151824c7587eafep+1)], [F(2048), F(0x1.a7d9a8edb47be8e4p+1)], 3585 [F(3), F(0x1.e8927964fd5fd08cp-2)], [F(5), F(0x1.65df657b04300868p-1)], 3586 [F(7), F(0x1.b0b0b0b78cc3f296p-1)], [F(15), F(0x1.2d145116c16ff856p+0)], 3587 [F(17), F(0x1.3afeb354b7d9731ap+0)], [F(31), F(0x1.7dc9e145867e62eap+0)], 3588 [F(33), F(0x1.84bd545e4baeddp+0)], [F(63), F(0x1.cca1950e4511e192p+0)], 3589 [F(65), F(0x1.d01b16f9433cf7b8p+0)], [F(-0), -F.infinity], [F(0), -F.infinity], 3590 [F(0.1), F(-0x1p+0)], [F(0.25), F(-0x1.34413509f79fef32p-1)], 3591 [F(0.75), F(-0x1.ffbfc2bbc7803758p-4)], [F(0.875), F(-0x1.db11ed766abf432ep-5)], 3592 [F(10), F(0x1p+0)], [F(100), F(0x1p+1)], [F(10000), F(0x1p+2)], 3593 ]; 3594 testLog10(vals); 3595 }} 3596 { 3597 float[2][16] vals = [ 3598 [float.nan, float.nan],[-float.nan, float.nan], 3599 [float.infinity, float.infinity], [-float.infinity, float.nan], 3600 [float.min_normal, -0x1.2f703p+5f], [-float.min_normal, float.nan], 3601 [float.max, 0x1.344136p+5f], [-float.max, float.nan], 3602 [float.min_normal / 2, -0x1.31d8b2p+5f], [-float.min_normal / 2, float.nan], 3603 [float.max / 2, 0x1.31d8b2p+5f], [-float.max / 2, float.nan], 3604 [float.min_normal / 3, -0x1.334156p+5f], [-float.min_normal / 3, float.nan], 3605 [float.max / 3, 0x1.30701p+5f], [-float.max / 3, float.nan], 3606 ]; 3607 testLog10(vals); 3608 } 3609 { 3610 double[2][16] vals = [ 3611 [double.nan, double.nan],[-double.nan, double.nan], 3612 [double.infinity, double.infinity], [-double.infinity, double.nan], 3613 [double.min_normal, -0x1.33a7146f72a42p+8], [-double.min_normal, double.nan], 3614 [double.max, 0x1.34413509f79ffp+8], [-double.max, double.nan], 3615 [double.min_normal / 2, -0x1.33f424bcb522p+8], [-double.min_normal / 2, double.nan], 3616 [double.max / 2, 0x1.33f424bcb522p+8], [-double.max / 2, double.nan], 3617 [double.min_normal / 3, -0x1.3421390dcbe37p+8], [-double.min_normal / 3, double.nan], 3618 [double.max / 3, 0x1.33c7106b9e609p+8], [-double.max / 3, double.nan], 3619 ]; 3620 testLog10(vals); 3621 } 3622 alias F = floatTraits!real; 3623 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple) 3624 {{ 3625 real[2][16] vals = [ 3626 [real.nan, real.nan],[-real.nan, real.nan], 3627 [real.infinity, real.infinity], [-real.infinity, real.nan], 3628 [real.min_normal, -0x1.343793004f503232p+12L], [-real.min_normal, real.nan], 3629 [real.max, 0x1.34413509f79fef32p+12L], [-real.max, real.nan], 3630 [real.min_normal / 2, -0x1.343c6405237810b2p+12L], [-real.min_normal / 2, real.nan], 3631 [real.max / 2, 0x1.343c6405237810b2p+12L], [-real.max / 2, real.nan], 3632 [real.min_normal / 3, -0x1.343f354a34e427bp+12L], [-real.min_normal / 3, real.nan], 3633 [real.max / 3, 0x1.343992c0120bf9b2p+12L], [-real.max / 3, real.nan], 3634 ]; 3635 testLog10(vals); 3636 }} 3637 } 3638 3639 /** 3640 * Calculates the natural logarithm of 1 + x. 3641 * 3642 * For very small x, log1p(x) will be more accurate than 3643 * log(1 + x). 3644 * 3645 * $(TABLE_SV 3646 * $(TR $(TH x) $(TH log1p(x)) $(TH divide by 0?) $(TH invalid?)) 3647 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) $(TD no)) 3648 * $(TR $(TD -1.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) 3649 * $(TR $(TD $(LT)-1.0) $(TD -$(NAN)) $(TD no) $(TD yes)) 3650 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) 3651 * ) 3652 */ 3653 pragma(inline, true) 3654 real log1p(real x) @safe pure nothrow @nogc 3655 { 3656 version (INLINE_YL2X) 3657 { 3658 // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5, 3659 // ie if -0.29 <= x <= 0.414 3660 import std.math.constants : LN2; 3661 return (core.math.fabs(x) <= 0.25) ? core.math.yl2xp1(x, LN2) : core.math.yl2x(x+1, LN2); 3662 } 3663 else 3664 return log1pImpl(x); 3665 } 3666 3667 /// ditto 3668 pragma(inline, true) 3669 double log1p(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log1p(cast(real) x) : log1pImpl(x); } 3670 3671 /// ditto 3672 pragma(inline, true) 3673 float log1p(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log1p(cast(real) x) : log1pImpl(x); } 3674 3675 // @@@DEPRECATED_[2.112.0]@@@ 3676 deprecated("`std.math.exponential.log1p` called with argument types `(int)` matches both " 3677 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.") 3678 real log1p(int x) @safe pure nothrow @nogc { return log1p(cast(real) x); } 3679 // @@@DEPRECATED_[2.112.0]@@@ 3680 deprecated("`std.math.exponential.log1p` called with argument types `(uint)` matches both " 3681 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.") 3682 real log1p(uint x) @safe pure nothrow @nogc { return log1p(cast(real) x); } 3683 // @@@DEPRECATED_[2.112.0]@@@ 3684 deprecated("`std.math.exponential.log1p` called with argument types `(long)` matches both " 3685 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.") 3686 real log1p(long x) @safe pure nothrow @nogc { return log1p(cast(real) x); } 3687 // @@@DEPRECATED_[2.112.0]@@@ 3688 deprecated("`std.math.exponential.log1p` called with argument types `(ulong)` matches both " 3689 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.") 3690 real log1p(ulong x) @safe pure nothrow @nogc { return log1p(cast(real) x); } 3691 3692 /// 3693 @safe pure unittest 3694 { 3695 import std.math.traits : isIdentical, isNaN; 3696 import std.math.operations : feqrel; 3697 3698 assert(isIdentical(log1p(0.0), 0.0)); 3699 assert(log1p(1.0).feqrel(0.69314) > 16); 3700 3701 assert(log1p(-1.0) == -real.infinity); 3702 assert(isNaN(log1p(-2.0))); 3703 assert(log1p(real.nan) is real.nan); 3704 assert(log1p(-real.nan) is -real.nan); 3705 assert(log1p(real.infinity) == real.infinity); 3706 } 3707 3708 private T log1pImpl(T)(T x) @safe pure nothrow @nogc 3709 { 3710 import std.math.traits : isNaN, isInfinity, signbit; 3711 import std.math.algebraic : poly; 3712 import std.math.constants : SQRT1_2, SQRT2; 3713 import std.math : floatTraits, RealFormat; 3714 3715 // Special cases. 3716 if (isNaN(x) || x == 0.0) 3717 return x; 3718 if (isInfinity(x) && !signbit(x)) 3719 return x; 3720 if (x == -1.0) 3721 return -T.infinity; 3722 if (x < -1.0) 3723 return T.nan; 3724 3725 alias F = floatTraits!T; 3726 static if (F.realFormat == RealFormat.ieeeSingle || 3727 F.realFormat == RealFormat.ieeeDouble) 3728 { 3729 // When the input is within the range 1/sqrt(2) <= x+1 <= sqrt(2), compute 3730 // log1p inline. Forwarding to log() would otherwise result in inaccuracies. 3731 const T xp1 = x + 1.0; 3732 if (xp1 >= SQRT1_2 && xp1 <= SQRT2) 3733 { 3734 alias coeffs = LogCoeffs!T; 3735 3736 T px = poly(x, coeffs.logp1P); 3737 T qx = poly(x, coeffs.logp1Q); 3738 const T xx = x * x; 3739 qx = x + ((cast(T) -0.5) * xx + x * (xx * px / qx)); 3740 return qx; 3741 } 3742 } 3743 3744 return logImpl!(T, true)(x); 3745 } 3746 3747 @safe @nogc nothrow unittest 3748 { 3749 import std.math : floatTraits, RealFormat; 3750 import std.meta : AliasSeq; 3751 3752 static void testLog1p(T)(T[2][] vals) 3753 { 3754 import std.math.operations : isClose; 3755 import std.math.traits : isNaN; 3756 foreach (ref pair; vals) 3757 { 3758 if (isNaN(pair[1])) 3759 assert(isNaN(log1p(pair[0]))); 3760 else 3761 assert(isClose(log1p(pair[0]), pair[1])); 3762 } 3763 } 3764 static foreach (F; AliasSeq!(float, double, real)) 3765 {{ 3766 scope F[2][] vals = [ 3767 [F(1), F(0x1.62e42fefa39ef358p-1)], [F(2), F(0x1.193ea7aad030a976p+0)], 3768 [F(4), F(0x1.9c041f7ed8d336bp+0)], [F(8), F(0x1.193ea7aad030a976p+1)], 3769 [F(16), F(0x1.6aa6bc1fa7f79cfp+1)], [F(32), F(0x1.bf8d8f4d5b8d1038p+1)], 3770 [F(64), F(0x1.0b292939429755ap+2)], [F(128), F(0x1.37072a9b5b6cb31p+2)], 3771 [F(256), F(0x1.63241004e9010ad8p+2)], [F(512), F(0x1.8f60adf041bde2a8p+2)], 3772 [F(1024), F(0x1.bbad39ebe1cc08b6p+2)], [F(2048), F(0x1.e801c1698ba4395cp+2)], 3773 [F(3), F(0x1.62e42fefa39ef358p+0)], [F(5), F(0x1.cab0bfa2a2002322p+0)], 3774 [F(7), F(0x1.0a2b23f3bab73682p+1)], [F(15), F(0x1.62e42fefa39ef358p+1)], 3775 [F(17), F(0x1.71f7b3a6b918664cp+1)], [F(31), F(0x1.bb9d3beb8c86b02ep+1)], 3776 [F(33), F(0x1.c35fc81b90df59c6p+1)], [F(63), F(0x1.0a2b23f3bab73682p+2)], 3777 [F(65), F(0x1.0c234da4a23a6686p+2)], [F(-0), F(-0x0p+0)], [F(0), F(0x0p+0)], 3778 [F(0.1), F(0x1.8663f793c46c69cp-4)], [F(0.25), F(0x1.c8ff7c79a9a21ac4p-3)], 3779 [F(0.75), F(0x1.1e85f5e7040d03ep-1)], [F(0.875), F(0x1.41d8fe84672ae646p-1)], 3780 [F(10), F(0x1.32ee3b77f374bb7cp+1)], [F(100), F(0x1.275e2271bba30be4p+2)], 3781 [F(10000), F(0x1.26bbed6fbd84182ep+3)], 3782 ]; 3783 testLog1p(vals); 3784 }} 3785 { 3786 float[2][16] vals = [ 3787 [float.nan, float.nan],[-float.nan, float.nan], 3788 [float.infinity, float.infinity], [-float.infinity, float.nan], 3789 [float.min_normal, 0x1p-126f], [-float.min_normal, -0x1p-126f], 3790 [float.max, 0x1.62e43p+6f], [-float.max, float.nan], 3791 [float.min_normal / 2, 0x0.8p-126f], [-float.min_normal / 2, -0x0.8p-126f], 3792 [float.max / 2, 0x1.601e68p+6f], [-float.max / 2, float.nan], 3793 [float.min_normal / 3, 0x0.555556p-126f], [-float.min_normal / 3, -0x0.555556p-126f], 3794 [float.max / 3, 0x1.5e7f36p+6f], [-float.max / 3, float.nan], 3795 ]; 3796 testLog1p(vals); 3797 } 3798 { 3799 double[2][16] vals = [ 3800 [double.nan, double.nan],[-double.nan, double.nan], 3801 [double.infinity, double.infinity], [-double.infinity, double.nan], 3802 [double.min_normal, 0x1p-1022], [-double.min_normal, -0x1p-1022], 3803 [double.max, 0x1.62e42fefa39efp+9], [-double.max, double.nan], 3804 [double.min_normal / 2, 0x0.8p-1022], [-double.min_normal / 2, -0x0.8p-1022], 3805 [double.max / 2, 0x1.628b76e3a7b61p+9], [-double.max / 2, double.nan], 3806 [double.min_normal / 3, 0x0.5555555555555p-1022], [-double.min_normal / 3, -0x0.5555555555555p-1022], 3807 [double.max / 3, 0x1.6257909bce36ep+9], [-double.max / 3, double.nan], 3808 ]; 3809 testLog1p(vals); 3810 } 3811 alias F = floatTraits!real; 3812 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple) 3813 {{ 3814 real[2][16] vals = [ 3815 [real.nan, real.nan],[-real.nan, real.nan], 3816 [real.infinity, real.infinity], [-real.infinity, real.nan], 3817 [real.min_normal, 0x1p-16382L], [-real.min_normal, -0x1p-16382L], 3818 [real.max, 0x1.62e42fefa39ef358p+13L], [-real.max, real.nan], 3819 [real.min_normal / 2, 0x0.8p-16382L], [-real.min_normal / 2, -0x0.8p-16382L], 3820 [real.max / 2, 0x1.62dea45ee3e064dcp+13L], [-real.max / 2, real.nan], 3821 [real.min_normal / 3, 0x0.5555555555555556p-16382L], [-real.min_normal / 3, -0x0.5555555555555556p-16382L], 3822 [real.max / 3, 0x1.62db65fa664871d2p+13L], [-real.max / 3, real.nan], 3823 ]; 3824 testLog1p(vals); 3825 }} 3826 } 3827 3828 /*************************************** 3829 * Calculates the base-2 logarithm of x: 3830 * $(SUB log, 2)x 3831 * 3832 * $(TABLE_SV 3833 * $(TR $(TH x) $(TH log2(x)) $(TH divide by 0?) $(TH invalid?)) 3834 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no) ) 3835 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes) ) 3836 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no) ) 3837 * ) 3838 */ 3839 pragma(inline, true) 3840 real log2(real x) @safe pure nothrow @nogc 3841 { 3842 version (INLINE_YL2X) 3843 return core.math.yl2x(x, 1.0L); 3844 else 3845 return log2Impl(x); 3846 } 3847 3848 /// ditto 3849 pragma(inline, true) 3850 double log2(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log2(cast(real) x) : log2Impl(x); } 3851 3852 /// ditto 3853 pragma(inline, true) 3854 float log2(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log2(cast(real) x) : log2Impl(x); } 3855 3856 // @@@DEPRECATED_[2.112.0]@@@ 3857 deprecated("`std.math.exponential.log2` called with argument types `(int)` matches both " 3858 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.") 3859 real log2(int x) @safe pure nothrow @nogc { return log2(cast(real) x); } 3860 // @@@DEPRECATED_[2.112.0]@@@ 3861 deprecated("`std.math.exponential.log2` called with argument types `(uint)` matches both " 3862 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.") 3863 real log2(uint x) @safe pure nothrow @nogc { return log2(cast(real) x); } 3864 // @@@DEPRECATED_[2.112.0]@@@ 3865 deprecated("`std.math.exponential.log2` called with argument types `(long)` matches both " 3866 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.") 3867 real log2(long x) @safe pure nothrow @nogc { return log2(cast(real) x); } 3868 // @@@DEPRECATED_[2.112.0]@@@ 3869 deprecated("`std.math.exponential.log2` called with argument types `(ulong)` matches both " 3870 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.") 3871 real log2(ulong x) @safe pure nothrow @nogc { return log2(cast(real) x); } 3872 3873 /// 3874 @safe unittest 3875 { 3876 import std.math.operations : isClose; 3877 3878 assert(isClose(log2(1024.0L), 10)); 3879 } 3880 3881 @safe @nogc nothrow unittest 3882 { 3883 import std.math.operations : isClose; 3884 3885 // check if values are equal to 19 decimal digits of precision 3886 assert(isClose(log2(1024.0L), 10, 1e-18)); 3887 } 3888 3889 private T log2Impl(T)(T x) @safe pure nothrow @nogc 3890 { 3891 import std.math.traits : isNaN, isInfinity, signbit; 3892 import std.math.constants : SQRT1_2, LOG2E; 3893 import std.math.algebraic : poly; 3894 import std.math : floatTraits, RealFormat; 3895 3896 alias coeffs = LogCoeffs!T; 3897 alias F = floatTraits!T; 3898 3899 // Special cases are the same as for log. 3900 if (isNaN(x)) 3901 return x; 3902 if (isInfinity(x) && !signbit(x)) 3903 return x; 3904 if (x == 0.0) 3905 return -T.infinity; 3906 if (x < 0.0) 3907 return T.nan; 3908 3909 // Separate mantissa from exponent. 3910 // Note, frexp is used so that denormal numbers will be handled properly. 3911 T y, z; 3912 int exp; 3913 3914 x = frexp(x, exp); 3915 3916 static if (F.realFormat == RealFormat.ieeeDouble || 3917 F.realFormat == RealFormat.ieeeExtended || 3918 F.realFormat == RealFormat.ieeeExtended53 || 3919 F.realFormat == RealFormat.ieeeQuadruple) 3920 { 3921 // Logarithm using log(x) = z + z^^3 R(z) / S(z), 3922 // where z = 2(x - 1)/(x + 1) 3923 if ((exp > 2) || (exp < -2)) 3924 { 3925 if (x < SQRT1_2) 3926 { // 2(2x - 1)/(2x + 1) 3927 exp -= 1; 3928 z = x - 0.5; 3929 y = 0.5 * z + 0.5; 3930 } 3931 else 3932 { // 2(x - 1)/(x + 1) 3933 z = x - 0.5; 3934 z -= 0.5; 3935 y = 0.5 * x + 0.5; 3936 } 3937 x = z / y; 3938 z = x * x; 3939 y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS)); 3940 goto Ldone; 3941 } 3942 } 3943 3944 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x) 3945 if (x < SQRT1_2) 3946 { 3947 exp -= 1; 3948 x = 2.0 * x - 1.0; 3949 } 3950 else 3951 x = x - 1.0; 3952 3953 z = x * x; 3954 static if (F.realFormat == RealFormat.ieeeSingle) 3955 y = x * (z * poly(x, coeffs.log2P)); 3956 else 3957 y = x * (z * poly(x, coeffs.log2P) / poly(x, coeffs.log2Q)); 3958 y = y - 0.5 * z; 3959 3960 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2). 3961 // This sequence of operations is critical and it may be horribly 3962 // defeated by some compiler optimizers. 3963 Ldone: 3964 z = y * (LOG2E - 1.0); 3965 z += x * (LOG2E - 1.0); 3966 z += y; 3967 z += x; 3968 z += exp; 3969 3970 return z; 3971 } 3972 3973 @safe @nogc nothrow unittest 3974 { 3975 import std.math : floatTraits, RealFormat; 3976 import std.meta : AliasSeq; 3977 3978 static void testLog2(T)(T[2][] vals) 3979 { 3980 import std.math.operations : isClose; 3981 import std.math.traits : isNaN; 3982 foreach (ref pair; vals) 3983 { 3984 if (isNaN(pair[1])) 3985 assert(isNaN(log2(pair[0]))); 3986 else 3987 assert(isClose(log2(pair[0]), pair[1])); 3988 } 3989 } 3990 static foreach (F; AliasSeq!(float, double, real)) 3991 {{ 3992 scope F[2][] vals = [ 3993 [F(1), F(0x0p+0)], [F(2), F(0x1p+0)], 3994 [F(4), F(0x1p+1)], [F(8), F(0x1.8p+1)], 3995 [F(16), F(0x1p+2)], [F(32), F(0x1.4p+2)], 3996 [F(64), F(0x1.8p+2)], [F(128), F(0x1.cp+2)], 3997 [F(256), F(0x1p+3)], [F(512), F(0x1.2p+3)], 3998 [F(1024), F(0x1.4p+3)], [F(2048), F(0x1.6p+3)], 3999 [F(3), F(0x1.95c01a39fbd687ap+0)], [F(5), F(0x1.2934f0979a3715fcp+1)], 4000 [F(7), F(0x1.675767f54042cd9ap+1)], [F(15), F(0x1.f414fdb4982259ccp+1)], 4001 [F(17), F(0x1.0598fdbeb244c5ap+2)], [F(31), F(0x1.3d118d66c4d4e554p+2)], 4002 [F(33), F(0x1.42d75a6eb1dfb0e6p+2)], [F(63), F(0x1.7e8bc1179e0caa9cp+2)], 4003 [F(65), F(0x1.816e79685c2d2298p+2)], [F(-0), -F.infinity], [F(0), -F.infinity], 4004 [F(0.1), F(-0x1.a934f0979a3715fcp+1)], [F(0.25), F(-0x1p+1)], 4005 [F(0.75), F(-0x1.a8ff971810a5e182p-2)], [F(0.875), F(-0x1.8a8980abfbd32668p-3)], 4006 [F(10), F(0x1.a934f0979a3715fcp+1)], [F(100), F(0x1.a934f0979a3715fcp+2)], 4007 [F(10000), F(0x1.a934f0979a3715fcp+3)], 4008 ]; 4009 testLog2(vals); 4010 }} 4011 { 4012 float[2][16] vals = [ 4013 [float.nan, float.nan],[-float.nan, float.nan], 4014 [float.infinity, float.infinity], [-float.infinity, float.nan], 4015 [float.min_normal, -0x1.f8p+6f], [-float.min_normal, float.nan], 4016 [float.max, 0x1p+7f], [-float.max, float.nan], 4017 [float.min_normal / 2, -0x1.fcp+6f], [-float.min_normal / 2, float.nan], 4018 [float.max / 2, 0x1.fcp+6f], [-float.max / 2, float.nan], 4019 [float.min_normal / 3, -0x1.fe57p+6f], [-float.min_normal / 3, float.nan], 4020 [float.max / 3, 0x1.f9a9p+6f], [-float.max / 3, float.nan], 4021 ]; 4022 testLog2(vals); 4023 } 4024 { 4025 double[2][16] vals = [ 4026 [double.nan, double.nan],[-double.nan, double.nan], 4027 [double.infinity, double.infinity], [-double.infinity, double.nan], 4028 [double.min_normal, -0x1.ffp+9], [-double.min_normal, double.nan], 4029 [double.max, 0x1p+10], [-double.max, double.nan], 4030 [double.min_normal / 2, -0x1.ff8p+9], [-double.min_normal / 2, double.nan], 4031 [double.max / 2, 0x1.ff8p+9], [-double.max / 2, double.nan], 4032 [double.min_normal / 3, -0x1.ffcae00d1cfdfp+9], [-double.min_normal / 3, double.nan], 4033 [double.max / 3, 0x1.ff351ff2e3021p+9], [-double.max / 3, double.nan], 4034 ]; 4035 testLog2(vals); 4036 } 4037 alias F = floatTraits!real; 4038 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple) 4039 {{ 4040 real[2][16] vals = [ 4041 [real.nan, real.nan],[-real.nan, real.nan], 4042 [real.infinity, real.infinity], [-real.infinity, real.nan], 4043 [real.min_normal, -0x1.fffp+13L], [-real.min_normal, real.nan], 4044 [real.max, 0x1p+14L], [-real.max, real.nan], 4045 [real.min_normal / 2, -0x1.fff8p+13L], [-real.min_normal / 2, real.nan], 4046 [real.max / 2, 0x1.fff8p+13L], [-real.max / 2, real.nan], 4047 [real.min_normal / 3, -0x1.fffcae00d1cfdeb4p+13L], [-real.min_normal / 3, real.nan], 4048 [real.max / 3, 0x1.fff351ff2e30214cp+13L], [-real.max / 3, real.nan], 4049 ]; 4050 testLog2(vals); 4051 }} 4052 } 4053 4054 /***************************************** 4055 * Extracts the exponent of x as a signed integral value. 4056 * 4057 * If x is subnormal, it is treated as if it were normalized. 4058 * For a positive, finite x: 4059 * 4060 * 1 $(LT)= $(I x) * FLT_RADIX$(SUPERSCRIPT -logb(x)) $(LT) FLT_RADIX 4061 * 4062 * $(TABLE_SV 4063 * $(TR $(TH x) $(TH logb(x)) $(TH divide by 0?) ) 4064 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no)) 4065 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) ) 4066 * ) 4067 */ 4068 pragma(inline, true) 4069 real logb(real x) @trusted pure nothrow @nogc 4070 { 4071 version (InlineAsm_X87_MSVC) 4072 return logbAsm(x); 4073 else 4074 return logbImpl(x); 4075 } 4076 4077 /// ditto 4078 pragma(inline, true) 4079 double logb(double x) @trusted pure nothrow @nogc { return logbImpl(x); } 4080 4081 /// ditto 4082 pragma(inline, true) 4083 float logb(float x) @trusted pure nothrow @nogc { return logbImpl(x); } 4084 4085 /// 4086 @safe @nogc nothrow unittest 4087 { 4088 assert(logb(1.0) == 0); 4089 assert(logb(100.0) == 6); 4090 4091 assert(logb(0.0) == -real.infinity); 4092 assert(logb(real.infinity) == real.infinity); 4093 assert(logb(-real.infinity) == real.infinity); 4094 } 4095 4096 @safe @nogc nothrow unittest 4097 { 4098 import std.meta : AliasSeq; 4099 import std.typecons : Tuple; 4100 import std.math.traits : isNaN; 4101 static foreach (F; AliasSeq!(float, double, real)) 4102 {{ 4103 alias T = Tuple!(F, F); 4104 T[17] vals = // x, logb(x) 4105 [ 4106 T(1.0 , 0 ), 4107 T(100.0 , 6 ), 4108 T(0.0 , -F.infinity), 4109 T(-0.0 , -F.infinity), 4110 T(1024 , 10 ), 4111 T(-2000 , 10 ), 4112 T(0x0.1p-127 , -131 ), 4113 T(0x0.01p-127 , -135 ), 4114 T(0x0.011p-127 , -135 ), 4115 T(F.nan , F.nan ), 4116 T(-F.nan , F.nan ), 4117 T(F.infinity , F.infinity ), 4118 T(-F.infinity , F.infinity ), 4119 T(F.min_normal , F.min_exp-1), 4120 T(-F.min_normal, F.min_exp-1), 4121 T(F.max , F.max_exp-1), 4122 T(-F.max , F.max_exp-1), 4123 ]; 4124 4125 foreach (elem; vals) 4126 { 4127 if (isNaN(elem[1])) 4128 assert(isNaN(logb(elem[1]))); 4129 else 4130 assert(logb(elem[0]) == elem[1]); 4131 } 4132 }} 4133 } 4134 4135 version (InlineAsm_X87_MSVC) 4136 private T logbAsm(T)(T x) @trusted pure nothrow @nogc 4137 { 4138 version (X86_64) 4139 { 4140 asm pure nothrow @nogc 4141 { 4142 naked ; 4143 fld real ptr [RCX] ; 4144 fxtract ; 4145 fstp ST(0) ; 4146 ret ; 4147 } 4148 } 4149 else 4150 { 4151 asm pure nothrow @nogc 4152 { 4153 fld x ; 4154 fxtract ; 4155 fstp ST(0) ; 4156 } 4157 } 4158 } 4159 4160 private T logbImpl(T)(T x) @trusted pure nothrow @nogc 4161 { 4162 import std.math.traits : isFinite; 4163 4164 // Handle special cases. 4165 if (!isFinite(x)) 4166 return x * x; 4167 if (x == 0) 4168 return -1 / (x * x); 4169 4170 return ilogb(x); 4171 } 4172 4173 @safe @nogc nothrow unittest 4174 { 4175 import std.math : floatTraits, RealFormat; 4176 import std.meta : AliasSeq; 4177 4178 static void testLogb(T)(T[2][] vals) 4179 { 4180 import std.math.operations : isClose; 4181 import std.math.traits : isNaN; 4182 foreach (ref pair; vals) 4183 { 4184 if (isNaN(pair[1])) 4185 assert(isNaN(logb(pair[0]))); 4186 else 4187 assert(isClose(logb(pair[0]), pair[1])); 4188 } 4189 } 4190 static foreach (F; AliasSeq!(float, double, real)) 4191 {{ 4192 scope F[2][] vals = [ 4193 [F(1), F(0x0p+0)], [F(2), F(0x1p+0)], 4194 [F(4), F(0x1p+1)], [F(8), F(0x1.8p+1)], 4195 [F(16), F(0x1p+2)], [F(32), F(0x1.4p+2)], 4196 [F(64), F(0x1.8p+2)], [F(128), F(0x1.cp+2)], 4197 [F(256), F(0x1p+3)], [F(512), F(0x1.2p+3)], 4198 [F(1024), F(0x1.4p+3)], [F(2048), F(0x1.6p+3)], 4199 [F(3), F(0x1p+0)], [F(5), F(0x1p+1)], 4200 [F(7), F(0x1p+1)], [F(15), F(0x1.8p+1)], 4201 [F(17), F(0x1p+2)], [F(31), F(0x1p+2)], 4202 [F(33), F(0x1.4p+2)], [F(63), F(0x1.4p+2)], 4203 [F(65), F(0x1.8p+2)], [F(-0), -F.infinity], [F(0), -F.infinity], 4204 [F(0.1), F(-0x1p+2)], [F(0.25), F(-0x1p+1)], 4205 [F(0.75), F(-0x1p+0)], [F(0.875), F(-0x1p+0)], 4206 [F(10), F(0x1.8p+1)], [F(100), F(0x1.8p+2)], 4207 [F(10000), F(0x1.ap+3)], 4208 ]; 4209 testLogb(vals); 4210 }} 4211 { 4212 float[2][16] vals = [ 4213 [float.nan, float.nan],[-float.nan, float.nan], 4214 [float.infinity, float.infinity], [-float.infinity, float.infinity], 4215 [float.min_normal, -0x1.f8p+6f], [-float.min_normal, -0x1.f8p+6f], 4216 [float.max, 0x1.fcp+6f], [-float.max, 0x1.fcp+6f], 4217 [float.min_normal / 2, -0x1.fcp+6f], [-float.min_normal / 2, -0x1.fcp+6f], 4218 [float.max / 2, 0x1.f8p+6f], [-float.max / 2, 0x1.f8p+6f], 4219 [float.min_normal / 3, -0x1p+7f], [-float.min_normal / 3, -0x1p+7f], 4220 [float.max / 3, 0x1.f8p+6f], [-float.max / 3, 0x1.f8p+6f], 4221 ]; 4222 testLogb(vals); 4223 } 4224 { 4225 double[2][16] vals = [ 4226 [double.nan, double.nan],[-double.nan, double.nan], 4227 [double.infinity, double.infinity], [-double.infinity, double.infinity], 4228 [double.min_normal, -0x1.ffp+9], [-double.min_normal, -0x1.ffp+9], 4229 [double.max, 0x1.ff8p+9], [-double.max, 0x1.ff8p+9], 4230 [double.min_normal / 2, -0x1.ff8p+9], [-double.min_normal / 2, -0x1.ff8p+9], 4231 [double.max / 2, 0x1.ffp+9], [-double.max / 2, 0x1.ffp+9], 4232 [double.min_normal / 3, -0x1p+10], [-double.min_normal / 3, -0x1p+10], 4233 [double.max / 3, 0x1.ffp+9], [-double.max / 3, 0x1.ffp+9], 4234 ]; 4235 testLogb(vals); 4236 } 4237 alias F = floatTraits!real; 4238 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple) 4239 {{ 4240 real[2][16] vals = [ 4241 [real.nan, real.nan],[-real.nan, real.nan], 4242 [real.infinity, real.infinity], [-real.infinity, real.infinity], 4243 [real.min_normal, -0x1.fffp+13L], [-real.min_normal, -0x1.fffp+13L], 4244 [real.max, 0x1.fff8p+13L], [-real.max, 0x1.fff8p+13L], 4245 [real.min_normal / 2, -0x1.fff8p+13L], [-real.min_normal / 2, -0x1.fff8p+13L], 4246 [real.max / 2, 0x1.fffp+13L], [-real.max / 2, 0x1.fffp+13L], 4247 [real.min_normal / 3, -0x1p+14L], [-real.min_normal / 3, -0x1p+14L], 4248 [real.max / 3, 0x1.fffp+13L], [-real.max / 3, 0x1.fffp+13L], 4249 ]; 4250 testLogb(vals); 4251 }} 4252 } 4253 4254 /************************************* 4255 * Efficiently calculates x * 2$(SUPERSCRIPT n). 4256 * 4257 * scalbn handles underflow and overflow in 4258 * the same fashion as the basic arithmetic operators. 4259 * 4260 * $(TABLE_SV 4261 * $(TR $(TH x) $(TH scalb(x))) 4262 * $(TR $(TD $(PLUSMNINF)) $(TD $(PLUSMNINF)) ) 4263 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) ) 4264 * ) 4265 */ 4266 pragma(inline, true) 4267 real scalbn(real x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); } 4268 4269 /// ditto 4270 pragma(inline, true) 4271 double scalbn(double x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); } 4272 4273 /// ditto 4274 pragma(inline, true) 4275 float scalbn(float x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); } 4276 4277 /// 4278 @safe pure nothrow @nogc unittest 4279 { 4280 assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L); 4281 assert(scalbn(-real.infinity, 5) == -real.infinity); 4282 assert(scalbn(2.0,10) == 2048.0); 4283 assert(scalbn(2048.0f,-10) == 2.0f); 4284 } 4285 4286 pragma(inline, true) 4287 private F _scalbn(F)(F x, int n) 4288 { 4289 import std.math.traits : isInfinity; 4290 4291 if (__ctfe) 4292 { 4293 // Handle special cases. 4294 if (x == F(0.0) || isInfinity(x)) 4295 return x; 4296 } 4297 return core.math.ldexp(x, n); 4298 } 4299 4300 @safe pure nothrow @nogc unittest 4301 { 4302 // CTFE-able test 4303 static assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L); 4304 static assert(scalbn(-real.infinity, 5) == -real.infinity); 4305 // Test with large exponent delta n where the result is in bounds but 2.0L ^^ n is not. 4306 enum initialExponent = real.min_exp + 2, resultExponent = real.max_exp - 2; 4307 enum int n = resultExponent - initialExponent; 4308 enum real x = 0x1.2345678abcdefp0L * (2.0L ^^ initialExponent); 4309 enum staticResult = scalbn(x, n); 4310 static assert(staticResult == 0x1.2345678abcdefp0L * (2.0L ^^ resultExponent)); 4311 assert(scalbn(x, n) == staticResult); 4312 } 4313