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