1 // Written in the D programming language. 2 3 /** 4 This is a submodule of $(MREF std, math). 5 6 It contains several functions for work with floating point numbers. 7 8 Copyright: Copyright The D Language Foundation 2000 - 2011. 9 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). 10 Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston, 11 Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger 12 Source: $(PHOBOSSRC std/math/operations.d) 13 14 Macros: 15 TABLE_SV = <table border="1" cellpadding="4" cellspacing="0"> 16 <caption>Special Values</caption> 17 $0</table> 18 SVH = $(TR $(TH $1) $(TH $2)) 19 SV = $(TR $(TD $1) $(TD $2)) 20 NAN = $(RED NAN) 21 PLUSMN = ± 22 INFIN = ∞ 23 LT = < 24 GT = > 25 */ 26 27 module std.math.operations; 28 29 import std.traits : CommonType, isFloatingPoint, isIntegral, Unqual; 30 31 // Functions for NaN payloads 32 /* 33 * A 'payload' can be stored in the significand of a $(NAN). One bit is required 34 * to distinguish between a quiet and a signalling $(NAN). This leaves 22 bits 35 * of payload for a float; 51 bits for a double; 62 bits for an 80-bit real; 36 * and 111 bits for a 128-bit quad. 37 */ 38 /** 39 * Create a quiet $(NAN), storing an integer inside the payload. 40 * 41 * For floats, the largest possible payload is 0x3F_FFFF. 42 * For doubles, it is 0x3_FFFF_FFFF_FFFF. 43 * For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF. 44 */ 45 real NaN(ulong payload) @trusted pure nothrow @nogc 46 { 47 import std.math : floatTraits, RealFormat; 48 49 alias F = floatTraits!(real); 50 static if (F.realFormat == RealFormat.ieeeExtended || 51 F.realFormat == RealFormat.ieeeExtended53) 52 { 53 // real80 (in x86 real format, the implied bit is actually 54 // not implied but a real bit which is stored in the real) 55 ulong v = 3; // implied bit = 1, quiet bit = 1 56 } 57 else 58 { 59 ulong v = 1; // no implied bit. quiet bit = 1 60 } 61 if (__ctfe) 62 { 63 v = 1; // We use a double in CTFE. 64 assert(payload >>> 51 == 0, 65 "Cannot set more than 51 bits of NaN payload in CTFE."); 66 } 67 68 69 ulong a = payload; 70 71 // 22 Float bits 72 ulong w = a & 0x3F_FFFF; 73 a -= w; 74 75 v <<=22; 76 v |= w; 77 a >>=22; 78 79 // 29 Double bits 80 v <<=29; 81 w = a & 0xFFF_FFFF; 82 v |= w; 83 a -= w; 84 a >>=29; 85 86 if (__ctfe) 87 { 88 v |= 0x7FF0_0000_0000_0000; 89 return *cast(double*) &v; 90 } 91 else static if (F.realFormat == RealFormat.ieeeDouble) 92 { 93 v |= 0x7FF0_0000_0000_0000; 94 real x; 95 * cast(ulong *)(&x) = v; 96 return x; 97 } 98 else 99 { 100 v <<=11; 101 a &= 0x7FF; 102 v |= a; 103 real x = real.nan; 104 105 // Extended real bits 106 static if (F.realFormat == RealFormat.ieeeQuadruple) 107 { 108 v <<= 1; // there's no implicit bit 109 110 version (LittleEndian) 111 { 112 *cast(ulong*)(6+cast(ubyte*)(&x)) = v; 113 } 114 else 115 { 116 *cast(ulong*)(2+cast(ubyte*)(&x)) = v; 117 } 118 } 119 else 120 { 121 *cast(ulong *)(&x) = v; 122 } 123 return x; 124 } 125 } 126 127 /// 128 @safe @nogc pure nothrow unittest 129 { 130 import std.math.traits : isNaN; 131 132 real a = NaN(1_000_000); 133 assert(isNaN(a)); 134 assert(getNaNPayload(a) == 1_000_000); 135 } 136 137 @system pure nothrow @nogc unittest // not @safe because taking address of local. 138 { 139 import std.math : floatTraits, RealFormat; 140 141 static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) 142 { 143 auto x = NaN(1); 144 auto xl = *cast(ulong*)&x; 145 assert(xl & 0x8_0000_0000_0000UL); //non-signaling bit, bit 52 146 assert((xl & 0x7FF0_0000_0000_0000UL) == 0x7FF0_0000_0000_0000UL); //all exp bits set 147 } 148 } 149 150 /** 151 * Extract an integral payload from a $(NAN). 152 * 153 * Returns: 154 * the integer payload as a ulong. 155 * 156 * For floats, the largest possible payload is 0x3F_FFFF. 157 * For doubles, it is 0x3_FFFF_FFFF_FFFF. 158 * For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF. 159 */ 160 ulong getNaNPayload(real x) @trusted pure nothrow @nogc 161 { 162 import std.math : floatTraits, RealFormat; 163 164 // assert(isNaN(x)); 165 alias F = floatTraits!(real); 166 ulong m = void; 167 if (__ctfe) 168 { 169 double y = x; 170 m = *cast(ulong*) &y; 171 // Make it look like an 80-bit significand. 172 // Skip exponent, and quiet bit 173 m &= 0x0007_FFFF_FFFF_FFFF; 174 m <<= 11; 175 } 176 else static if (F.realFormat == RealFormat.ieeeDouble) 177 { 178 m = *cast(ulong*)(&x); 179 // Make it look like an 80-bit significand. 180 // Skip exponent, and quiet bit 181 m &= 0x0007_FFFF_FFFF_FFFF; 182 m <<= 11; 183 } 184 else static if (F.realFormat == RealFormat.ieeeQuadruple) 185 { 186 version (LittleEndian) 187 { 188 m = *cast(ulong*)(6+cast(ubyte*)(&x)); 189 } 190 else 191 { 192 m = *cast(ulong*)(2+cast(ubyte*)(&x)); 193 } 194 195 m >>= 1; // there's no implicit bit 196 } 197 else 198 { 199 m = *cast(ulong*)(&x); 200 } 201 202 // ignore implicit bit and quiet bit 203 204 const ulong f = m & 0x3FFF_FF00_0000_0000L; 205 206 ulong w = f >>> 40; 207 w |= (m & 0x00FF_FFFF_F800L) << (22 - 11); 208 w |= (m & 0x7FF) << 51; 209 return w; 210 } 211 212 /// 213 @safe @nogc pure nothrow unittest 214 { 215 import std.math.traits : isNaN; 216 217 real a = NaN(1_000_000); 218 assert(isNaN(a)); 219 assert(getNaNPayload(a) == 1_000_000); 220 } 221 222 @safe @nogc pure nothrow unittest 223 { 224 import std.math.traits : isIdentical, isNaN; 225 226 enum real a = NaN(1_000_000); 227 static assert(isNaN(a)); 228 static assert(getNaNPayload(a) == 1_000_000); 229 real b = NaN(1_000_000); 230 assert(isIdentical(b, a)); 231 // The CTFE version of getNaNPayload relies on it being impossible 232 // for a CTFE-constructed NaN to have more than 51 bits of payload. 233 enum nanNaN = NaN(getNaNPayload(real.nan)); 234 assert(isIdentical(real.nan, nanNaN)); 235 static if (real.init != real.init) 236 { 237 enum initNaN = NaN(getNaNPayload(real.init)); 238 assert(isIdentical(real.init, initNaN)); 239 } 240 } 241 242 debug(UnitTest) 243 { 244 @safe pure nothrow @nogc unittest 245 { 246 real nan4 = NaN(0x789_ABCD_EF12_3456); 247 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended 248 || floatTraits!(real).realFormat == RealFormat.ieeeQuadruple) 249 { 250 assert(getNaNPayload(nan4) == 0x789_ABCD_EF12_3456); 251 } 252 else 253 { 254 assert(getNaNPayload(nan4) == 0x1_ABCD_EF12_3456); 255 } 256 double nan5 = nan4; 257 assert(getNaNPayload(nan5) == 0x1_ABCD_EF12_3456); 258 float nan6 = nan4; 259 assert(getNaNPayload(nan6) == 0x12_3456); 260 nan4 = NaN(0xFABCD); 261 assert(getNaNPayload(nan4) == 0xFABCD); 262 nan6 = nan4; 263 assert(getNaNPayload(nan6) == 0xFABCD); 264 nan5 = NaN(0x100_0000_0000_3456); 265 assert(getNaNPayload(nan5) == 0x0000_0000_3456); 266 } 267 } 268 269 /** 270 * Calculate the next largest floating point value after x. 271 * 272 * Return the least number greater than x that is representable as a real; 273 * thus, it gives the next point on the IEEE number line. 274 * 275 * $(TABLE_SV 276 * $(SVH x, nextUp(x) ) 277 * $(SV -$(INFIN), -real.max ) 278 * $(SV $(PLUSMN)0.0, real.min_normal*real.epsilon ) 279 * $(SV real.max, $(INFIN) ) 280 * $(SV $(INFIN), $(INFIN) ) 281 * $(SV $(NAN), $(NAN) ) 282 * ) 283 */ 284 real nextUp(real x) @trusted pure nothrow @nogc 285 { 286 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB; 287 288 alias F = floatTraits!(real); 289 static if (F.realFormat != RealFormat.ieeeDouble) 290 { 291 if (__ctfe) 292 { 293 if (x == -real.infinity) 294 return -real.max; 295 if (!(x < real.infinity)) // Infinity or NaN. 296 return x; 297 real delta; 298 // Start with a decent estimate of delta. 299 if (x <= 0x1.ffffffffffffep+1023 && x >= -double.max) 300 { 301 const double d = cast(double) x; 302 delta = (cast(real) nextUp(d) - cast(real) d) * 0x1p-11L; 303 while (x + (delta * 0x1p-100L) > x) 304 delta *= 0x1p-100L; 305 } 306 else 307 { 308 delta = 0x1p960L; 309 while (!(x + delta > x) && delta < real.max * 0x1p-100L) 310 delta *= 0x1p100L; 311 } 312 if (x + delta > x) 313 { 314 while (x + (delta / 2) > x) 315 delta /= 2; 316 } 317 else 318 { 319 do { delta += delta; } while (!(x + delta > x)); 320 } 321 if (x < 0 && x + delta == 0) 322 return -0.0L; 323 return x + delta; 324 } 325 } 326 static if (F.realFormat == RealFormat.ieeeDouble) 327 { 328 return nextUp(cast(double) x); 329 } 330 else static if (F.realFormat == RealFormat.ieeeQuadruple) 331 { 332 ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]; 333 if (e == F.EXPMASK) 334 { 335 // NaN or Infinity 336 if (x == -real.infinity) return -real.max; 337 return x; // +Inf and NaN are unchanged. 338 } 339 340 auto ps = cast(ulong *)&x; 341 if (ps[MANTISSA_MSB] & 0x8000_0000_0000_0000) 342 { 343 // Negative number 344 if (ps[MANTISSA_LSB] == 0 && ps[MANTISSA_MSB] == 0x8000_0000_0000_0000) 345 { 346 // it was negative zero, change to smallest subnormal 347 ps[MANTISSA_LSB] = 1; 348 ps[MANTISSA_MSB] = 0; 349 return x; 350 } 351 if (ps[MANTISSA_LSB] == 0) --ps[MANTISSA_MSB]; 352 --ps[MANTISSA_LSB]; 353 } 354 else 355 { 356 // Positive number 357 ++ps[MANTISSA_LSB]; 358 if (ps[MANTISSA_LSB] == 0) ++ps[MANTISSA_MSB]; 359 } 360 return x; 361 } 362 else static if (F.realFormat == RealFormat.ieeeExtended || 363 F.realFormat == RealFormat.ieeeExtended53) 364 { 365 // For 80-bit reals, the "implied bit" is a nuisance... 366 ushort *pe = cast(ushort *)&x; 367 ulong *ps = cast(ulong *)&x; 368 // EPSILON is 1 for 64-bit, and 2048 for 53-bit precision reals. 369 enum ulong EPSILON = 2UL ^^ (64 - real.mant_dig); 370 371 if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK) 372 { 373 // First, deal with NANs and infinity 374 if (x == -real.infinity) return -real.max; 375 return x; // +Inf and NaN are unchanged. 376 } 377 if (pe[F.EXPPOS_SHORT] & 0x8000) 378 { 379 // Negative number -- need to decrease the significand 380 *ps -= EPSILON; 381 // Need to mask with 0x7FFF... so subnormals are treated correctly. 382 if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF) 383 { 384 if (pe[F.EXPPOS_SHORT] == 0x8000) // it was negative zero 385 { 386 *ps = 1; 387 pe[F.EXPPOS_SHORT] = 0; // smallest subnormal. 388 return x; 389 } 390 391 --pe[F.EXPPOS_SHORT]; 392 393 if (pe[F.EXPPOS_SHORT] == 0x8000) 394 return x; // it's become a subnormal, implied bit stays low. 395 396 *ps = 0xFFFF_FFFF_FFFF_FFFF; // set the implied bit 397 return x; 398 } 399 return x; 400 } 401 else 402 { 403 // Positive number -- need to increase the significand. 404 // Works automatically for positive zero. 405 *ps += EPSILON; 406 if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0) 407 { 408 // change in exponent 409 ++pe[F.EXPPOS_SHORT]; 410 *ps = 0x8000_0000_0000_0000; // set the high bit 411 } 412 } 413 return x; 414 } 415 else // static if (F.realFormat == RealFormat.ibmExtended) 416 { 417 assert(0, "nextUp not implemented"); 418 } 419 } 420 421 /** ditto */ 422 double nextUp(double x) @trusted pure nothrow @nogc 423 { 424 ulong s = *cast(ulong *)&x; 425 426 if ((s & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000) 427 { 428 // First, deal with NANs and infinity 429 if (x == -x.infinity) return -x.max; 430 return x; // +INF and NAN are unchanged. 431 } 432 if (s & 0x8000_0000_0000_0000) // Negative number 433 { 434 if (s == 0x8000_0000_0000_0000) // it was negative zero 435 { 436 s = 0x0000_0000_0000_0001; // change to smallest subnormal 437 return *cast(double*) &s; 438 } 439 --s; 440 } 441 else 442 { // Positive number 443 ++s; 444 } 445 return *cast(double*) &s; 446 } 447 448 /** ditto */ 449 float nextUp(float x) @trusted pure nothrow @nogc 450 { 451 uint s = *cast(uint *)&x; 452 453 if ((s & 0x7F80_0000) == 0x7F80_0000) 454 { 455 // First, deal with NANs and infinity 456 if (x == -x.infinity) return -x.max; 457 458 return x; // +INF and NAN are unchanged. 459 } 460 if (s & 0x8000_0000) // Negative number 461 { 462 if (s == 0x8000_0000) // it was negative zero 463 { 464 s = 0x0000_0001; // change to smallest subnormal 465 return *cast(float*) &s; 466 } 467 468 --s; 469 } 470 else 471 { 472 // Positive number 473 ++s; 474 } 475 return *cast(float*) &s; 476 } 477 478 /// 479 @safe @nogc pure nothrow unittest 480 { 481 assert(nextUp(1.0 - 1.0e-6).feqrel(0.999999) > 16); 482 assert(nextUp(1.0 - real.epsilon).feqrel(1.0) > 16); 483 } 484 485 /** 486 * Calculate the next smallest floating point value before x. 487 * 488 * Return the greatest number less than x that is representable as a real; 489 * thus, it gives the previous point on the IEEE number line. 490 * 491 * $(TABLE_SV 492 * $(SVH x, nextDown(x) ) 493 * $(SV $(INFIN), real.max ) 494 * $(SV $(PLUSMN)0.0, -real.min_normal*real.epsilon ) 495 * $(SV -real.max, -$(INFIN) ) 496 * $(SV -$(INFIN), -$(INFIN) ) 497 * $(SV $(NAN), $(NAN) ) 498 * ) 499 */ 500 real nextDown(real x) @safe pure nothrow @nogc 501 { 502 return -nextUp(-x); 503 } 504 505 /** ditto */ 506 double nextDown(double x) @safe pure nothrow @nogc 507 { 508 return -nextUp(-x); 509 } 510 511 /** ditto */ 512 float nextDown(float x) @safe pure nothrow @nogc 513 { 514 return -nextUp(-x); 515 } 516 517 /// 518 @safe pure nothrow @nogc unittest 519 { 520 assert( nextDown(1.0 + real.epsilon) == 1.0); 521 } 522 523 @safe pure nothrow @nogc unittest 524 { 525 import std.math : floatTraits, RealFormat; 526 import std.math.traits : isIdentical; 527 528 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended || 529 floatTraits!(real).realFormat == RealFormat.ieeeDouble || 530 floatTraits!(real).realFormat == RealFormat.ieeeExtended53 || 531 floatTraits!(real).realFormat == RealFormat.ieeeQuadruple) 532 { 533 // Tests for reals 534 assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC))); 535 //static assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC))); 536 // negative numbers 537 assert( nextUp(-real.infinity) == -real.max ); 538 assert( nextUp(-1.0L-real.epsilon) == -1.0 ); 539 assert( nextUp(-2.0L) == -2.0 + real.epsilon); 540 static assert( nextUp(-real.infinity) == -real.max ); 541 static assert( nextUp(-1.0L-real.epsilon) == -1.0 ); 542 static assert( nextUp(-2.0L) == -2.0 + real.epsilon); 543 // subnormals and zero 544 assert( nextUp(-real.min_normal) == -real.min_normal*(1-real.epsilon) ); 545 assert( nextUp(-real.min_normal*(1-real.epsilon)) == -real.min_normal*(1-2*real.epsilon) ); 546 assert( isIdentical(-0.0L, nextUp(-real.min_normal*real.epsilon)) ); 547 assert( nextUp(-0.0L) == real.min_normal*real.epsilon ); 548 assert( nextUp(0.0L) == real.min_normal*real.epsilon ); 549 assert( nextUp(real.min_normal*(1-real.epsilon)) == real.min_normal ); 550 assert( nextUp(real.min_normal) == real.min_normal*(1+real.epsilon) ); 551 static assert( nextUp(-real.min_normal) == -real.min_normal*(1-real.epsilon) ); 552 static assert( nextUp(-real.min_normal*(1-real.epsilon)) == -real.min_normal*(1-2*real.epsilon) ); 553 static assert( -0.0L is nextUp(-real.min_normal*real.epsilon) ); 554 static assert( nextUp(-0.0L) == real.min_normal*real.epsilon ); 555 static assert( nextUp(0.0L) == real.min_normal*real.epsilon ); 556 static assert( nextUp(real.min_normal*(1-real.epsilon)) == real.min_normal ); 557 static assert( nextUp(real.min_normal) == real.min_normal*(1+real.epsilon) ); 558 // positive numbers 559 assert( nextUp(1.0L) == 1.0 + real.epsilon ); 560 assert( nextUp(2.0L-real.epsilon) == 2.0 ); 561 assert( nextUp(real.max) == real.infinity ); 562 assert( nextUp(real.infinity)==real.infinity ); 563 static assert( nextUp(1.0L) == 1.0 + real.epsilon ); 564 static assert( nextUp(2.0L-real.epsilon) == 2.0 ); 565 static assert( nextUp(real.max) == real.infinity ); 566 static assert( nextUp(real.infinity)==real.infinity ); 567 // ctfe near double.max boundary 568 static assert(nextUp(nextDown(cast(real) double.max)) == cast(real) double.max); 569 } 570 571 double n = NaN(0xABC); 572 assert(isIdentical(nextUp(n), n)); 573 // negative numbers 574 assert( nextUp(-double.infinity) == -double.max ); 575 assert( nextUp(-1-double.epsilon) == -1.0 ); 576 assert( nextUp(-2.0) == -2.0 + double.epsilon); 577 // subnormals and zero 578 579 assert( nextUp(-double.min_normal) == -double.min_normal*(1-double.epsilon) ); 580 assert( nextUp(-double.min_normal*(1-double.epsilon)) == -double.min_normal*(1-2*double.epsilon) ); 581 assert( isIdentical(-0.0, nextUp(-double.min_normal*double.epsilon)) ); 582 assert( nextUp(0.0) == double.min_normal*double.epsilon ); 583 assert( nextUp(-0.0) == double.min_normal*double.epsilon ); 584 assert( nextUp(double.min_normal*(1-double.epsilon)) == double.min_normal ); 585 assert( nextUp(double.min_normal) == double.min_normal*(1+double.epsilon) ); 586 // positive numbers 587 assert( nextUp(1.0) == 1.0 + double.epsilon ); 588 assert( nextUp(2.0-double.epsilon) == 2.0 ); 589 assert( nextUp(double.max) == double.infinity ); 590 591 float fn = NaN(0xABC); 592 assert(isIdentical(nextUp(fn), fn)); 593 float f = -float.min_normal*(1-float.epsilon); 594 float f1 = -float.min_normal; 595 assert( nextUp(f1) == f); 596 f = 1.0f+float.epsilon; 597 f1 = 1.0f; 598 assert( nextUp(f1) == f ); 599 f1 = -0.0f; 600 assert( nextUp(f1) == float.min_normal*float.epsilon); 601 assert( nextUp(float.infinity)==float.infinity ); 602 603 assert(nextDown(1.0L+real.epsilon)==1.0); 604 assert(nextDown(1.0+double.epsilon)==1.0); 605 f = 1.0f+float.epsilon; 606 assert(nextDown(f)==1.0); 607 assert(nextafter(1.0+real.epsilon, -real.infinity)==1.0); 608 609 // CTFE 610 611 enum double ctfe_n = NaN(0xABC); 612 //static assert(isIdentical(nextUp(ctfe_n), ctfe_n)); // FIXME: https://issues.dlang.org/show_bug.cgi?id=20197 613 static assert(nextUp(double.nan) is double.nan); 614 // negative numbers 615 static assert( nextUp(-double.infinity) == -double.max ); 616 static assert( nextUp(-1-double.epsilon) == -1.0 ); 617 static assert( nextUp(-2.0) == -2.0 + double.epsilon); 618 // subnormals and zero 619 620 static assert( nextUp(-double.min_normal) == -double.min_normal*(1-double.epsilon) ); 621 static assert( nextUp(-double.min_normal*(1-double.epsilon)) == -double.min_normal*(1-2*double.epsilon) ); 622 static assert( -0.0 is nextUp(-double.min_normal*double.epsilon) ); 623 static assert( nextUp(0.0) == double.min_normal*double.epsilon ); 624 static assert( nextUp(-0.0) == double.min_normal*double.epsilon ); 625 static assert( nextUp(double.min_normal*(1-double.epsilon)) == double.min_normal ); 626 static assert( nextUp(double.min_normal) == double.min_normal*(1+double.epsilon) ); 627 // positive numbers 628 static assert( nextUp(1.0) == 1.0 + double.epsilon ); 629 static assert( nextUp(2.0-double.epsilon) == 2.0 ); 630 static assert( nextUp(double.max) == double.infinity ); 631 632 enum float ctfe_fn = NaN(0xABC); 633 //static assert(isIdentical(nextUp(ctfe_fn), ctfe_fn)); // FIXME: https://issues.dlang.org/show_bug.cgi?id=20197 634 static assert(nextUp(float.nan) is float.nan); 635 static assert(nextUp(-float.min_normal) == -float.min_normal*(1-float.epsilon)); 636 static assert(nextUp(1.0f) == 1.0f+float.epsilon); 637 static assert(nextUp(-0.0f) == float.min_normal*float.epsilon); 638 static assert(nextUp(float.infinity)==float.infinity); 639 static assert(nextDown(1.0L+real.epsilon)==1.0); 640 static assert(nextDown(1.0+double.epsilon)==1.0); 641 static assert(nextDown(1.0f+float.epsilon)==1.0); 642 static assert(nextafter(1.0+real.epsilon, -real.infinity)==1.0); 643 } 644 645 646 647 /****************************************** 648 * Calculates the next representable value after x in the direction of y. 649 * 650 * If y > x, the result will be the next largest floating-point value; 651 * if y < x, the result will be the next smallest value. 652 * If x == y, the result is y. 653 * If x or y is a NaN, the result is a NaN. 654 * 655 * Remarks: 656 * This function is not generally very useful; it's almost always better to use 657 * the faster functions nextUp() or nextDown() instead. 658 * 659 * The FE_INEXACT and FE_OVERFLOW exceptions will be raised if x is finite and 660 * the function result is infinite. The FE_INEXACT and FE_UNDERFLOW 661 * exceptions will be raised if the function value is subnormal, and x is 662 * not equal to y. 663 */ 664 T nextafter(T)(const T x, const T y) @safe pure nothrow @nogc 665 { 666 import std.math.traits : isNaN; 667 668 if (x == y || isNaN(y)) 669 { 670 return y; 671 } 672 673 if (isNaN(x)) 674 { 675 return x; 676 } 677 678 return ((y>x) ? nextUp(x) : nextDown(x)); 679 } 680 681 /// 682 @safe pure nothrow @nogc unittest 683 { 684 import std.math.traits : isNaN; 685 686 float a = 1; 687 assert(is(typeof(nextafter(a, a)) == float)); 688 assert(nextafter(a, a.infinity) > a); 689 assert(isNaN(nextafter(a, a.nan))); 690 assert(isNaN(nextafter(a.nan, a))); 691 692 double b = 2; 693 assert(is(typeof(nextafter(b, b)) == double)); 694 assert(nextafter(b, b.infinity) > b); 695 assert(isNaN(nextafter(b, b.nan))); 696 assert(isNaN(nextafter(b.nan, b))); 697 698 real c = 3; 699 assert(is(typeof(nextafter(c, c)) == real)); 700 assert(nextafter(c, c.infinity) > c); 701 assert(isNaN(nextafter(c, c.nan))); 702 assert(isNaN(nextafter(c.nan, c))); 703 } 704 705 @safe pure nothrow @nogc unittest 706 { 707 import std.math.traits : isNaN, signbit; 708 709 // CTFE 710 enum float a = 1; 711 static assert(is(typeof(nextafter(a, a)) == float)); 712 static assert(nextafter(a, a.infinity) > a); 713 static assert(isNaN(nextafter(a, a.nan))); 714 static assert(isNaN(nextafter(a.nan, a))); 715 716 enum double b = 2; 717 static assert(is(typeof(nextafter(b, b)) == double)); 718 static assert(nextafter(b, b.infinity) > b); 719 static assert(isNaN(nextafter(b, b.nan))); 720 static assert(isNaN(nextafter(b.nan, b))); 721 722 enum real c = 3; 723 static assert(is(typeof(nextafter(c, c)) == real)); 724 static assert(nextafter(c, c.infinity) > c); 725 static assert(isNaN(nextafter(c, c.nan))); 726 static assert(isNaN(nextafter(c.nan, c))); 727 728 enum real negZero = nextafter(+0.0L, -0.0L); 729 static assert(negZero == -0.0L); 730 static assert(signbit(negZero)); 731 732 static assert(nextafter(c, c) == c); 733 } 734 735 //real nexttoward(real x, real y) { return core.stdc.math.nexttowardl(x, y); } 736 737 /** 738 * Returns the positive difference between x and y. 739 * 740 * Equivalent to `fmax(x-y, 0)`. 741 * 742 * Returns: 743 * $(TABLE_SV 744 * $(TR $(TH x, y) $(TH fdim(x, y))) 745 * $(TR $(TD x $(GT) y) $(TD x - y)) 746 * $(TR $(TD x $(LT)= y) $(TD +0.0)) 747 * ) 748 */ 749 real fdim(real x, real y) @safe pure nothrow @nogc 750 { 751 return (x < y) ? +0.0 : x - y; 752 } 753 754 /// 755 @safe pure nothrow @nogc unittest 756 { 757 import std.math.traits : isNaN; 758 759 assert(fdim(2.0, 0.0) == 2.0); 760 assert(fdim(-2.0, 0.0) == 0.0); 761 assert(fdim(real.infinity, 2.0) == real.infinity); 762 assert(isNaN(fdim(real.nan, 2.0))); 763 assert(isNaN(fdim(2.0, real.nan))); 764 assert(isNaN(fdim(real.nan, real.nan))); 765 } 766 767 /** 768 * Returns the larger of `x` and `y`. 769 * 770 * If one of the arguments is a `NaN`, the other is returned. 771 * 772 * See_Also: $(REF max, std,algorithm,comparison) is faster because it does not perform the `isNaN` test. 773 */ 774 F fmax(F)(const F x, const F y) @safe pure nothrow @nogc 775 if (__traits(isFloating, F)) 776 { 777 import std.math.traits : isNaN; 778 779 // Do the more predictable test first. Generates 0 branches with ldc and 1 branch with gdc. 780 // See https://godbolt.org/z/erxrW9 781 if (isNaN(x)) return y; 782 return y > x ? y : x; 783 } 784 785 /// 786 @safe pure nothrow @nogc unittest 787 { 788 import std.meta : AliasSeq; 789 static foreach (F; AliasSeq!(float, double, real)) 790 { 791 assert(fmax(F(0.0), F(2.0)) == 2.0); 792 assert(fmax(F(-2.0), 0.0) == F(0.0)); 793 assert(fmax(F.infinity, F(2.0)) == F.infinity); 794 assert(fmax(F.nan, F(2.0)) == F(2.0)); 795 assert(fmax(F(2.0), F.nan) == F(2.0)); 796 } 797 } 798 799 /** 800 * Returns the smaller of `x` and `y`. 801 * 802 * If one of the arguments is a `NaN`, the other is returned. 803 * 804 * See_Also: $(REF min, std,algorithm,comparison) is faster because it does not perform the `isNaN` test. 805 */ 806 F fmin(F)(const F x, const F y) @safe pure nothrow @nogc 807 if (__traits(isFloating, F)) 808 { 809 import std.math.traits : isNaN; 810 811 // Do the more predictable test first. Generates 0 branches with ldc and 1 branch with gdc. 812 // See https://godbolt.org/z/erxrW9 813 if (isNaN(x)) return y; 814 return y < x ? y : x; 815 } 816 817 /// 818 @safe pure nothrow @nogc unittest 819 { 820 import std.meta : AliasSeq; 821 static foreach (F; AliasSeq!(float, double, real)) 822 { 823 assert(fmin(F(0.0), F(2.0)) == 0.0); 824 assert(fmin(F(-2.0), F(0.0)) == -2.0); 825 assert(fmin(F.infinity, F(2.0)) == 2.0); 826 assert(fmin(F.nan, F(2.0)) == 2.0); 827 assert(fmin(F(2.0), F.nan) == 2.0); 828 } 829 } 830 831 /************************************** 832 * Returns (x * y) + z, rounding only once according to the 833 * current rounding mode. 834 * 835 * BUGS: Not currently implemented - rounds twice. 836 */ 837 pragma(inline, true) 838 real fma(real x, real y, real z) @safe pure nothrow @nogc { return (x * y) + z; } 839 840 /// 841 @safe pure nothrow @nogc unittest 842 { 843 assert(fma(0.0, 2.0, 2.0) == 2.0); 844 assert(fma(2.0, 2.0, 2.0) == 6.0); 845 assert(fma(real.infinity, 2.0, 2.0) == real.infinity); 846 assert(fma(real.nan, 2.0, 2.0) is real.nan); 847 assert(fma(2.0, 2.0, real.nan) is real.nan); 848 } 849 850 /************************************** 851 * To what precision is x equal to y? 852 * 853 * Returns: the number of mantissa bits which are equal in x and y. 854 * eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision. 855 * 856 * $(TABLE_SV 857 * $(TR $(TH x) $(TH y) $(TH feqrel(x, y))) 858 * $(TR $(TD x) $(TD x) $(TD real.mant_dig)) 859 * $(TR $(TD x) $(TD $(GT)= 2*x) $(TD 0)) 860 * $(TR $(TD x) $(TD $(LT)= x/2) $(TD 0)) 861 * $(TR $(TD $(NAN)) $(TD any) $(TD 0)) 862 * $(TR $(TD any) $(TD $(NAN)) $(TD 0)) 863 * ) 864 */ 865 int feqrel(X)(const X x, const X y) @trusted pure nothrow @nogc 866 if (isFloatingPoint!(X)) 867 { 868 import std.math : floatTraits, RealFormat; 869 import core.math : fabs; 870 871 /* Public Domain. Author: Don Clugston, 18 Aug 2005. 872 */ 873 alias F = floatTraits!(X); 874 static if (F.realFormat == RealFormat.ieeeSingle 875 || F.realFormat == RealFormat.ieeeDouble 876 || F.realFormat == RealFormat.ieeeExtended 877 || F.realFormat == RealFormat.ieeeExtended53 878 || F.realFormat == RealFormat.ieeeQuadruple) 879 { 880 if (x == y) 881 return X.mant_dig; // ensure diff != 0, cope with INF. 882 883 Unqual!X diff = fabs(x - y); 884 885 ushort *pa = cast(ushort *)(&x); 886 ushort *pb = cast(ushort *)(&y); 887 ushort *pd = cast(ushort *)(&diff); 888 889 890 // The difference in abs(exponent) between x or y and abs(x-y) 891 // is equal to the number of significand bits of x which are 892 // equal to y. If negative, x and y have different exponents. 893 // If positive, x and y are equal to 'bitsdiff' bits. 894 // AND with 0x7FFF to form the absolute value. 895 // To avoid out-by-1 errors, we subtract 1 so it rounds down 896 // if the exponents were different. This means 'bitsdiff' is 897 // always 1 lower than we want, except that if bitsdiff == 0, 898 // they could have 0 or 1 bits in common. 899 900 int bitsdiff = ((( (pa[F.EXPPOS_SHORT] & F.EXPMASK) 901 + (pb[F.EXPPOS_SHORT] & F.EXPMASK) 902 - (1 << F.EXPSHIFT)) >> 1) 903 - (pd[F.EXPPOS_SHORT] & F.EXPMASK)) >> F.EXPSHIFT; 904 if ( (pd[F.EXPPOS_SHORT] & F.EXPMASK) == 0) 905 { // Difference is subnormal 906 // For subnormals, we need to add the number of zeros that 907 // lie at the start of diff's significand. 908 // We do this by multiplying by 2^^real.mant_dig 909 diff *= F.RECIP_EPSILON; 910 return bitsdiff + X.mant_dig - ((pd[F.EXPPOS_SHORT] & F.EXPMASK) >> F.EXPSHIFT); 911 } 912 913 if (bitsdiff > 0) 914 return bitsdiff + 1; // add the 1 we subtracted before 915 916 // Avoid out-by-1 errors when factor is almost 2. 917 if (bitsdiff == 0 918 && ((pa[F.EXPPOS_SHORT] ^ pb[F.EXPPOS_SHORT]) & F.EXPMASK) == 0) 919 { 920 return 1; 921 } else return 0; 922 } 923 else 924 { 925 static assert(false, "Not implemented for this architecture"); 926 } 927 } 928 929 /// 930 @safe pure unittest 931 { 932 assert(feqrel(2.0, 2.0) == 53); 933 assert(feqrel(2.0f, 2.0f) == 24); 934 assert(feqrel(2.0, double.nan) == 0); 935 936 // Test that numbers are within n digits of each 937 // other by testing if feqrel > n * log2(10) 938 939 // five digits 940 assert(feqrel(2.0, 2.00001) > 16); 941 // ten digits 942 assert(feqrel(2.0, 2.00000000001) > 33); 943 } 944 945 @safe pure nothrow @nogc unittest 946 { 947 void testFeqrel(F)() 948 { 949 // Exact equality 950 assert(feqrel(F.max, F.max) == F.mant_dig); 951 assert(feqrel!(F)(0.0, 0.0) == F.mant_dig); 952 assert(feqrel(F.infinity, F.infinity) == F.mant_dig); 953 954 // a few bits away from exact equality 955 F w=1; 956 for (int i = 1; i < F.mant_dig - 1; ++i) 957 { 958 assert(feqrel!(F)(1.0 + w * F.epsilon, 1.0) == F.mant_dig-i); 959 assert(feqrel!(F)(1.0 - w * F.epsilon, 1.0) == F.mant_dig-i); 960 assert(feqrel!(F)(1.0, 1 + (w-1) * F.epsilon) == F.mant_dig - i + 1); 961 w*=2; 962 } 963 964 assert(feqrel!(F)(1.5+F.epsilon, 1.5) == F.mant_dig-1); 965 assert(feqrel!(F)(1.5-F.epsilon, 1.5) == F.mant_dig-1); 966 assert(feqrel!(F)(1.5-F.epsilon, 1.5+F.epsilon) == F.mant_dig-2); 967 968 969 // Numbers that are close 970 assert(feqrel!(F)(0x1.Bp+84, 0x1.B8p+84) == 5); 971 assert(feqrel!(F)(0x1.8p+10, 0x1.Cp+10) == 2); 972 assert(feqrel!(F)(1.5 * (1 - F.epsilon), 1.0L) == 2); 973 assert(feqrel!(F)(1.5, 1.0) == 1); 974 assert(feqrel!(F)(2 * (1 - F.epsilon), 1.0L) == 1); 975 976 // Factors of 2 977 assert(feqrel(F.max, F.infinity) == 0); 978 assert(feqrel!(F)(2 * (1 - F.epsilon), 1.0L) == 1); 979 assert(feqrel!(F)(1.0, 2.0) == 0); 980 assert(feqrel!(F)(4.0, 1.0) == 0); 981 982 // Extreme inequality 983 assert(feqrel(F.nan, F.nan) == 0); 984 assert(feqrel!(F)(0.0L, -F.nan) == 0); 985 assert(feqrel(F.nan, F.infinity) == 0); 986 assert(feqrel(F.infinity, -F.infinity) == 0); 987 assert(feqrel(F.max, -F.max) == 0); 988 989 assert(feqrel(F.min_normal / 8, F.min_normal / 17) == 3); 990 991 const F Const = 2; 992 immutable F Immutable = 2; 993 auto Compiles = feqrel(Const, Immutable); 994 } 995 996 assert(feqrel(7.1824L, 7.1824L) == real.mant_dig); 997 998 testFeqrel!(real)(); 999 testFeqrel!(double)(); 1000 testFeqrel!(float)(); 1001 } 1002 1003 /** 1004 Computes whether a values is approximately equal to a reference value, 1005 admitting a maximum relative difference, and a maximum absolute difference. 1006 1007 Warning: 1008 This template is considered out-dated. It will be removed from 1009 Phobos in 2.106.0. Please use $(LREF isClose) instead. To achieve 1010 a similar behaviour to `approxEqual(a, b)` use 1011 `isClose(a, b, 1e-2, 1e-5)`. In case of comparing to 0.0, 1012 `isClose(a, b, 0.0, eps)` should be used, where `eps` 1013 represents the accepted deviation from 0.0." 1014 1015 Params: 1016 value = Value to compare. 1017 reference = Reference value. 1018 maxRelDiff = Maximum allowable difference relative to `reference`. 1019 Setting to 0.0 disables this check. Defaults to `1e-2`. 1020 maxAbsDiff = Maximum absolute difference. This is mainly usefull 1021 for comparing values to zero. Setting to 0.0 disables this check. 1022 Defaults to `1e-5`. 1023 1024 Returns: 1025 `true` if `value` is approximately equal to `reference` under 1026 either criterium. It is sufficient, when `value ` satisfies 1027 one of the two criteria. 1028 1029 If one item is a range, and the other is a single value, then 1030 the result is the logical and-ing of calling `approxEqual` on 1031 each element of the ranged item against the single item. If 1032 both items are ranges, then `approxEqual` returns `true` if 1033 and only if the ranges have the same number of elements and if 1034 `approxEqual` evaluates to `true` for each pair of elements. 1035 1036 See_Also: 1037 Use $(LREF feqrel) to get the number of equal bits in the mantissa. 1038 */ 1039 deprecated("approxEqual will be removed in 2.106.0. Please use isClose instead.") 1040 bool approxEqual(T, U, V)(T value, U reference, V maxRelDiff = 1e-2, V maxAbsDiff = 1e-5) 1041 { 1042 import core.math : fabs; 1043 import std.range.primitives : empty, front, isInputRange, popFront; 1044 static if (isInputRange!T) 1045 { 1046 static if (isInputRange!U) 1047 { 1048 // Two ranges 1049 for (;; value.popFront(), reference.popFront()) 1050 { 1051 if (value.empty) return reference.empty; 1052 if (reference.empty) return value.empty; 1053 if (!approxEqual(value.front, reference.front, maxRelDiff, maxAbsDiff)) 1054 return false; 1055 } 1056 } 1057 else static if (isIntegral!U) 1058 { 1059 // convert reference to real 1060 return approxEqual(value, real(reference), maxRelDiff, maxAbsDiff); 1061 } 1062 else 1063 { 1064 // value is range, reference is number 1065 for (; !value.empty; value.popFront()) 1066 { 1067 if (!approxEqual(value.front, reference, maxRelDiff, maxAbsDiff)) 1068 return false; 1069 } 1070 return true; 1071 } 1072 } 1073 else 1074 { 1075 static if (isInputRange!U) 1076 { 1077 // value is number, reference is range 1078 for (; !reference.empty; reference.popFront()) 1079 { 1080 if (!approxEqual(value, reference.front, maxRelDiff, maxAbsDiff)) 1081 return false; 1082 } 1083 return true; 1084 } 1085 else static if (isIntegral!T || isIntegral!U) 1086 { 1087 // convert both value and reference to real 1088 return approxEqual(real(value), real(reference), maxRelDiff, maxAbsDiff); 1089 } 1090 else 1091 { 1092 // two numbers 1093 //static assert(is(T : real) && is(U : real)); 1094 if (reference == 0) 1095 { 1096 return fabs(value) <= maxAbsDiff; 1097 } 1098 static if (is(typeof(value.infinity)) && is(typeof(reference.infinity))) 1099 { 1100 if (value == value.infinity && reference == reference.infinity || 1101 value == -value.infinity && reference == -reference.infinity) return true; 1102 } 1103 return fabs((value - reference) / reference) <= maxRelDiff 1104 || maxAbsDiff != 0 && fabs(value - reference) <= maxAbsDiff; 1105 } 1106 } 1107 } 1108 1109 deprecated @safe pure nothrow unittest 1110 { 1111 assert(approxEqual(1.0, 1.0099)); 1112 assert(!approxEqual(1.0, 1.011)); 1113 assert(approxEqual(0.00001, 0.0)); 1114 assert(!approxEqual(0.00002, 0.0)); 1115 1116 assert(approxEqual(3.0, [3, 3.01, 2.99])); // several reference values is strange 1117 assert(approxEqual([3, 3.01, 2.99], 3.0)); // better 1118 1119 float[] arr1 = [ 1.0, 2.0, 3.0 ]; 1120 double[] arr2 = [ 1.001, 1.999, 3 ]; 1121 assert(approxEqual(arr1, arr2)); 1122 } 1123 1124 deprecated @safe pure nothrow unittest 1125 { 1126 // relative comparison depends on reference, make sure proper 1127 // side is used when comparing range to single value. Based on 1128 // https://issues.dlang.org/show_bug.cgi?id=15763 1129 auto a = [2e-3 - 1e-5]; 1130 auto b = 2e-3 + 1e-5; 1131 assert(a[0].approxEqual(b)); 1132 assert(!b.approxEqual(a[0])); 1133 assert(a.approxEqual(b)); 1134 assert(!b.approxEqual(a)); 1135 } 1136 1137 deprecated @safe pure nothrow @nogc unittest 1138 { 1139 assert(!approxEqual(0.0,1e-15,1e-9,0.0)); 1140 assert(approxEqual(0.0,1e-15,1e-9,1e-9)); 1141 assert(!approxEqual(1.0,3.0,0.0,1.0)); 1142 1143 assert(approxEqual(1.00000000099,1.0,1e-9,0.0)); 1144 assert(!approxEqual(1.0000000011,1.0,1e-9,0.0)); 1145 } 1146 1147 deprecated @safe pure nothrow @nogc unittest 1148 { 1149 // maybe unintuitive behavior 1150 assert(approxEqual(1000.0,1010.0)); 1151 assert(approxEqual(9_090_000_000.0,9_000_000_000.0)); 1152 assert(approxEqual(0.0,1e30,1.0)); 1153 assert(approxEqual(0.00001,1e-30)); 1154 assert(!approxEqual(-1e-30,1e-30,1e-2,0.0)); 1155 } 1156 1157 deprecated @safe pure nothrow @nogc unittest 1158 { 1159 int a = 10; 1160 assert(approxEqual(10, a)); 1161 1162 assert(!approxEqual(3, 0)); 1163 assert(approxEqual(3, 3)); 1164 assert(approxEqual(3.0, 3)); 1165 assert(approxEqual(3, 3.0)); 1166 1167 assert(approxEqual(0.0,0.0)); 1168 assert(approxEqual(-0.0,0.0)); 1169 assert(approxEqual(0.0f,0.0)); 1170 } 1171 1172 deprecated @safe pure nothrow @nogc unittest 1173 { 1174 real num = real.infinity; 1175 assert(num == real.infinity); 1176 assert(approxEqual(num, real.infinity)); 1177 num = -real.infinity; 1178 assert(num == -real.infinity); 1179 assert(approxEqual(num, -real.infinity)); 1180 1181 assert(!approxEqual(1,real.nan)); 1182 assert(!approxEqual(real.nan,real.max)); 1183 assert(!approxEqual(real.nan,real.nan)); 1184 } 1185 1186 deprecated @safe pure nothrow unittest 1187 { 1188 assert(!approxEqual([1.0,2.0,3.0],[1.0,2.0])); 1189 assert(!approxEqual([1.0,2.0],[1.0,2.0,3.0])); 1190 1191 assert(approxEqual!(real[],real[])([],[])); 1192 assert(approxEqual(cast(real[])[],cast(real[])[])); 1193 } 1194 1195 1196 /** 1197 Computes whether two values are approximately equal, admitting a maximum 1198 relative difference, and a maximum absolute difference. 1199 1200 Params: 1201 lhs = First item to compare. 1202 rhs = Second item to compare. 1203 maxRelDiff = Maximum allowable relative difference. 1204 Setting to 0.0 disables this check. Default depends on the type of 1205 `lhs` and `rhs`: It is approximately half the number of decimal digits of 1206 precision of the smaller type. 1207 maxAbsDiff = Maximum absolute difference. This is mainly usefull 1208 for comparing values to zero. Setting to 0.0 disables this check. 1209 Defaults to `0.0`. 1210 1211 Returns: 1212 `true` if the two items are approximately equal under either criterium. 1213 It is sufficient, when `value ` satisfies one of the two criteria. 1214 1215 If one item is a range, and the other is a single value, then 1216 the result is the logical and-ing of calling `isClose` on 1217 each element of the ranged item against the single item. If 1218 both items are ranges, then `isClose` returns `true` if 1219 and only if the ranges have the same number of elements and if 1220 `isClose` evaluates to `true` for each pair of elements. 1221 1222 See_Also: 1223 Use $(LREF feqrel) to get the number of equal bits in the mantissa. 1224 */ 1225 bool isClose(T, U, V = CommonType!(FloatingPointBaseType!T,FloatingPointBaseType!U)) 1226 (T lhs, U rhs, V maxRelDiff = CommonDefaultFor!(T,U), V maxAbsDiff = 0.0) 1227 { 1228 import std.range.primitives : empty, front, isInputRange, popFront; 1229 import std.complex : Complex; 1230 static if (isInputRange!T) 1231 { 1232 static if (isInputRange!U) 1233 { 1234 // Two ranges 1235 for (;; lhs.popFront(), rhs.popFront()) 1236 { 1237 if (lhs.empty) return rhs.empty; 1238 if (rhs.empty) return lhs.empty; 1239 if (!isClose(lhs.front, rhs.front, maxRelDiff, maxAbsDiff)) 1240 return false; 1241 } 1242 } 1243 else 1244 { 1245 // lhs is range, rhs is number 1246 for (; !lhs.empty; lhs.popFront()) 1247 { 1248 if (!isClose(lhs.front, rhs, maxRelDiff, maxAbsDiff)) 1249 return false; 1250 } 1251 return true; 1252 } 1253 } 1254 else static if (isInputRange!U) 1255 { 1256 // lhs is number, rhs is range 1257 for (; !rhs.empty; rhs.popFront()) 1258 { 1259 if (!isClose(lhs, rhs.front, maxRelDiff, maxAbsDiff)) 1260 return false; 1261 } 1262 return true; 1263 } 1264 else static if (is(T TE == Complex!TE)) 1265 { 1266 static if (is(U UE == Complex!UE)) 1267 { 1268 // Two complex numbers 1269 return isClose(lhs.re, rhs.re, maxRelDiff, maxAbsDiff) 1270 && isClose(lhs.im, rhs.im, maxRelDiff, maxAbsDiff); 1271 } 1272 else 1273 { 1274 // lhs is complex, rhs is number 1275 return isClose(lhs.re, rhs, maxRelDiff, maxAbsDiff) 1276 && isClose(lhs.im, 0.0, maxRelDiff, maxAbsDiff); 1277 } 1278 } 1279 else static if (is(U UE == Complex!UE)) 1280 { 1281 // lhs is number, rhs is complex 1282 return isClose(lhs, rhs.re, maxRelDiff, maxAbsDiff) 1283 && isClose(0.0, rhs.im, maxRelDiff, maxAbsDiff); 1284 } 1285 else 1286 { 1287 // two numbers 1288 if (lhs == rhs) return true; 1289 1290 static if (is(typeof(lhs.infinity))) 1291 if (lhs == lhs.infinity || lhs == -lhs.infinity) 1292 return false; 1293 static if (is(typeof(rhs.infinity))) 1294 if (rhs == rhs.infinity || rhs == -rhs.infinity) 1295 return false; 1296 1297 import std.math.algebraic : abs; 1298 1299 auto diff = abs(lhs - rhs); 1300 1301 return diff <= maxRelDiff*abs(lhs) 1302 || diff <= maxRelDiff*abs(rhs) 1303 || diff <= maxAbsDiff; 1304 } 1305 } 1306 1307 /// 1308 @safe pure nothrow @nogc unittest 1309 { 1310 assert(isClose(1.0,0.999_999_999)); 1311 assert(isClose(0.001, 0.000_999_999_999)); 1312 assert(isClose(1_000_000_000.0,999_999_999.0)); 1313 1314 assert(isClose(17.123_456_789, 17.123_456_78)); 1315 assert(!isClose(17.123_456_789, 17.123_45)); 1316 1317 // use explicit 3rd parameter for less (or more) accuracy 1318 assert(isClose(17.123_456_789, 17.123_45, 1e-6)); 1319 assert(!isClose(17.123_456_789, 17.123_45, 1e-7)); 1320 1321 // use 4th parameter when comparing close to zero 1322 assert(!isClose(1e-100, 0.0)); 1323 assert(isClose(1e-100, 0.0, 0.0, 1e-90)); 1324 assert(!isClose(1e-10, -1e-10)); 1325 assert(isClose(1e-10, -1e-10, 0.0, 1e-9)); 1326 assert(!isClose(1e-300, 1e-298)); 1327 assert(isClose(1e-300, 1e-298, 0.0, 1e-200)); 1328 1329 // different default limits for different floating point types 1330 assert(isClose(1.0f, 0.999_99f)); 1331 assert(!isClose(1.0, 0.999_99)); 1332 static if (real.sizeof > double.sizeof) 1333 assert(!isClose(1.0L, 0.999_999_999L)); 1334 } 1335 1336 /// 1337 @safe pure nothrow unittest 1338 { 1339 assert(isClose([1.0, 2.0, 3.0], [0.999_999_999, 2.000_000_001, 3.0])); 1340 assert(!isClose([1.0, 2.0], [0.999_999_999, 2.000_000_001, 3.0])); 1341 assert(!isClose([1.0, 2.0, 3.0], [0.999_999_999, 2.000_000_001])); 1342 1343 assert(isClose([2.0, 1.999_999_999, 2.000_000_001], 2.0)); 1344 assert(isClose(2.0, [2.0, 1.999_999_999, 2.000_000_001])); 1345 } 1346 1347 @safe pure nothrow unittest 1348 { 1349 assert(!isClose([1.0, 2.0, 3.0], [0.999_999_999, 3.0, 3.0])); 1350 assert(!isClose([2.0, 1.999_999, 2.000_000_001], 2.0)); 1351 assert(!isClose(2.0, [2.0, 1.999_999_999, 2.000_000_999])); 1352 } 1353 1354 @safe pure nothrow @nogc unittest 1355 { 1356 immutable a = 1.00001f; 1357 const b = 1.000019; 1358 assert(isClose(a,b)); 1359 1360 assert(isClose(1.00001f,1.000019f)); 1361 assert(isClose(1.00001f,1.000019)); 1362 assert(isClose(1.00001,1.000019f)); 1363 assert(!isClose(1.00001,1.000019)); 1364 1365 real a1 = 1e-300L; 1366 real a2 = a1.nextUp; 1367 assert(isClose(a1,a2)); 1368 } 1369 1370 @safe pure nothrow unittest 1371 { 1372 float[] arr1 = [ 1.0, 2.0, 3.0 ]; 1373 double[] arr2 = [ 1.00001, 1.99999, 3 ]; 1374 assert(isClose(arr1, arr2)); 1375 } 1376 1377 @safe pure nothrow @nogc unittest 1378 { 1379 assert(!isClose(1000.0,1010.0)); 1380 assert(!isClose(9_090_000_000.0,9_000_000_000.0)); 1381 assert(isClose(0.0,1e30,1.0)); 1382 assert(!isClose(0.00001,1e-30)); 1383 assert(!isClose(-1e-30,1e-30,1e-2,0.0)); 1384 } 1385 1386 @safe pure nothrow @nogc unittest 1387 { 1388 assert(!isClose(3, 0)); 1389 assert(isClose(3, 3)); 1390 assert(isClose(3.0, 3)); 1391 assert(isClose(3, 3.0)); 1392 1393 assert(isClose(0.0,0.0)); 1394 assert(isClose(-0.0,0.0)); 1395 assert(isClose(0.0f,0.0)); 1396 } 1397 1398 @safe pure nothrow @nogc unittest 1399 { 1400 real num = real.infinity; 1401 assert(num == real.infinity); 1402 assert(isClose(num, real.infinity)); 1403 num = -real.infinity; 1404 assert(num == -real.infinity); 1405 assert(isClose(num, -real.infinity)); 1406 1407 assert(!isClose(1,real.nan)); 1408 assert(!isClose(real.nan,real.max)); 1409 assert(!isClose(real.nan,real.nan)); 1410 1411 assert(!isClose(-double.infinity, 1)); 1412 } 1413 1414 @safe pure nothrow @nogc unittest 1415 { 1416 assert(isClose!(real[],real[],real)([],[])); 1417 assert(isClose(cast(real[])[],cast(real[])[])); 1418 } 1419 1420 @safe pure nothrow @nogc unittest 1421 { 1422 import std.conv : to; 1423 1424 float f = 31.79f; 1425 double d = 31.79; 1426 double f2d = f.to!double; 1427 1428 assert(isClose(f,f2d)); 1429 assert(!isClose(d,f2d)); 1430 } 1431 1432 @safe pure nothrow @nogc unittest 1433 { 1434 import std.conv : to; 1435 1436 double d = 31.79; 1437 float f = d.to!float; 1438 double f2d = f.to!double; 1439 1440 assert(isClose(f,f2d)); 1441 assert(!isClose(d,f2d)); 1442 assert(isClose(d,f2d,1e-4)); 1443 } 1444 1445 package(std.math) template CommonDefaultFor(T,U) 1446 { 1447 import std.algorithm.comparison : min; 1448 1449 alias baseT = FloatingPointBaseType!T; 1450 alias baseU = FloatingPointBaseType!U; 1451 1452 enum CommonType!(baseT, baseU) CommonDefaultFor = 10.0L ^^ -((min(baseT.dig, baseU.dig) + 1) / 2 + 1); 1453 } 1454 1455 private template FloatingPointBaseType(T) 1456 { 1457 import std.range.primitives : ElementType; 1458 static if (isFloatingPoint!T) 1459 { 1460 alias FloatingPointBaseType = Unqual!T; 1461 } 1462 else static if (isFloatingPoint!(ElementType!(Unqual!T))) 1463 { 1464 alias FloatingPointBaseType = Unqual!(ElementType!(Unqual!T)); 1465 } 1466 else 1467 { 1468 alias FloatingPointBaseType = real; 1469 } 1470 } 1471 1472 /*********************************** 1473 * Defines a total order on all floating-point numbers. 1474 * 1475 * The order is defined as follows: 1476 * $(UL 1477 * $(LI All numbers in [-$(INFIN), +$(INFIN)] are ordered 1478 * the same way as by built-in comparison, with the exception of 1479 * -0.0, which is less than +0.0;) 1480 * $(LI If the sign bit is set (that is, it's 'negative'), $(NAN) is less 1481 * than any number; if the sign bit is not set (it is 'positive'), 1482 * $(NAN) is greater than any number;) 1483 * $(LI $(NAN)s of the same sign are ordered by the payload ('negative' 1484 * ones - in reverse order).) 1485 * ) 1486 * 1487 * Returns: 1488 * negative value if `x` precedes `y` in the order specified above; 1489 * 0 if `x` and `y` are identical, and positive value otherwise. 1490 * 1491 * See_Also: 1492 * $(MYREF isIdentical) 1493 * Standards: Conforms to IEEE 754-2008 1494 */ 1495 int cmp(T)(const(T) x, const(T) y) @nogc @trusted pure nothrow 1496 if (isFloatingPoint!T) 1497 { 1498 import std.math : floatTraits, RealFormat; 1499 1500 alias F = floatTraits!T; 1501 1502 static if (F.realFormat == RealFormat.ieeeSingle 1503 || F.realFormat == RealFormat.ieeeDouble) 1504 { 1505 static if (T.sizeof == 4) 1506 alias UInt = uint; 1507 else 1508 alias UInt = ulong; 1509 1510 union Repainter 1511 { 1512 T number; 1513 UInt bits; 1514 } 1515 1516 enum msb = ~(UInt.max >>> 1); 1517 1518 import std.typecons : Tuple; 1519 Tuple!(Repainter, Repainter) vars = void; 1520 vars[0].number = x; 1521 vars[1].number = y; 1522 1523 foreach (ref var; vars) 1524 if (var.bits & msb) 1525 var.bits = ~var.bits; 1526 else 1527 var.bits |= msb; 1528 1529 if (vars[0].bits < vars[1].bits) 1530 return -1; 1531 else if (vars[0].bits > vars[1].bits) 1532 return 1; 1533 else 1534 return 0; 1535 } 1536 else static if (F.realFormat == RealFormat.ieeeExtended53 1537 || F.realFormat == RealFormat.ieeeExtended 1538 || F.realFormat == RealFormat.ieeeQuadruple) 1539 { 1540 static if (F.realFormat == RealFormat.ieeeQuadruple) 1541 alias RemT = ulong; 1542 else 1543 alias RemT = ushort; 1544 1545 struct Bits 1546 { 1547 ulong bulk; 1548 RemT rem; 1549 } 1550 1551 union Repainter 1552 { 1553 T number; 1554 Bits bits; 1555 ubyte[T.sizeof] bytes; 1556 } 1557 1558 import std.typecons : Tuple; 1559 Tuple!(Repainter, Repainter) vars = void; 1560 vars[0].number = x; 1561 vars[1].number = y; 1562 1563 foreach (ref var; vars) 1564 if (var.bytes[F.SIGNPOS_BYTE] & 0x80) 1565 { 1566 var.bits.bulk = ~var.bits.bulk; 1567 var.bits.rem = cast(typeof(var.bits.rem))(-1 - var.bits.rem); // ~var.bits.rem 1568 } 1569 else 1570 { 1571 var.bytes[F.SIGNPOS_BYTE] |= 0x80; 1572 } 1573 1574 version (LittleEndian) 1575 { 1576 if (vars[0].bits.rem < vars[1].bits.rem) 1577 return -1; 1578 else if (vars[0].bits.rem > vars[1].bits.rem) 1579 return 1; 1580 else if (vars[0].bits.bulk < vars[1].bits.bulk) 1581 return -1; 1582 else if (vars[0].bits.bulk > vars[1].bits.bulk) 1583 return 1; 1584 else 1585 return 0; 1586 } 1587 else 1588 { 1589 if (vars[0].bits.bulk < vars[1].bits.bulk) 1590 return -1; 1591 else if (vars[0].bits.bulk > vars[1].bits.bulk) 1592 return 1; 1593 else if (vars[0].bits.rem < vars[1].bits.rem) 1594 return -1; 1595 else if (vars[0].bits.rem > vars[1].bits.rem) 1596 return 1; 1597 else 1598 return 0; 1599 } 1600 } 1601 else 1602 { 1603 // IBM Extended doubledouble does not follow the general 1604 // sign-exponent-significand layout, so has to be handled generically 1605 1606 import std.math.traits : signbit, isNaN; 1607 1608 const int xSign = signbit(x), 1609 ySign = signbit(y); 1610 1611 if (xSign == 1 && ySign == 1) 1612 return cmp(-y, -x); 1613 else if (xSign == 1) 1614 return -1; 1615 else if (ySign == 1) 1616 return 1; 1617 else if (x < y) 1618 return -1; 1619 else if (x == y) 1620 return 0; 1621 else if (x > y) 1622 return 1; 1623 else if (isNaN(x) && !isNaN(y)) 1624 return 1; 1625 else if (isNaN(y) && !isNaN(x)) 1626 return -1; 1627 else if (getNaNPayload(x) < getNaNPayload(y)) 1628 return -1; 1629 else if (getNaNPayload(x) > getNaNPayload(y)) 1630 return 1; 1631 else 1632 return 0; 1633 } 1634 } 1635 1636 /// Most numbers are ordered naturally. 1637 @safe unittest 1638 { 1639 assert(cmp(-double.infinity, -double.max) < 0); 1640 assert(cmp(-double.max, -100.0) < 0); 1641 assert(cmp(-100.0, -0.5) < 0); 1642 assert(cmp(-0.5, 0.0) < 0); 1643 assert(cmp(0.0, 0.5) < 0); 1644 assert(cmp(0.5, 100.0) < 0); 1645 assert(cmp(100.0, double.max) < 0); 1646 assert(cmp(double.max, double.infinity) < 0); 1647 1648 assert(cmp(1.0, 1.0) == 0); 1649 } 1650 1651 /// Positive and negative zeroes are distinct. 1652 @safe unittest 1653 { 1654 assert(cmp(-0.0, +0.0) < 0); 1655 assert(cmp(+0.0, -0.0) > 0); 1656 } 1657 1658 /// Depending on the sign, $(NAN)s go to either end of the spectrum. 1659 @safe unittest 1660 { 1661 assert(cmp(-double.nan, -double.infinity) < 0); 1662 assert(cmp(double.infinity, double.nan) < 0); 1663 assert(cmp(-double.nan, double.nan) < 0); 1664 } 1665 1666 /// $(NAN)s of the same sign are ordered by the payload. 1667 @safe unittest 1668 { 1669 assert(cmp(NaN(10), NaN(20)) < 0); 1670 assert(cmp(-NaN(20), -NaN(10)) < 0); 1671 } 1672 1673 @safe unittest 1674 { 1675 import std.meta : AliasSeq; 1676 static foreach (T; AliasSeq!(float, double, real)) 1677 {{ 1678 T[] values = [-cast(T) NaN(20), -cast(T) NaN(10), -T.nan, -T.infinity, 1679 -T.max, -T.max / 2, T(-16.0), T(-1.0).nextDown, 1680 T(-1.0), T(-1.0).nextUp, 1681 T(-0.5), -T.min_normal, (-T.min_normal).nextUp, 1682 -2 * T.min_normal * T.epsilon, 1683 -T.min_normal * T.epsilon, 1684 T(-0.0), T(0.0), 1685 T.min_normal * T.epsilon, 1686 2 * T.min_normal * T.epsilon, 1687 T.min_normal.nextDown, T.min_normal, T(0.5), 1688 T(1.0).nextDown, T(1.0), 1689 T(1.0).nextUp, T(16.0), T.max / 2, T.max, 1690 T.infinity, T.nan, cast(T) NaN(10), cast(T) NaN(20)]; 1691 1692 foreach (i, x; values) 1693 { 1694 foreach (y; values[i + 1 .. $]) 1695 { 1696 assert(cmp(x, y) < 0); 1697 assert(cmp(y, x) > 0); 1698 } 1699 assert(cmp(x, x) == 0); 1700 } 1701 }} 1702 } 1703 1704 package(std): // not yet public 1705 1706 struct FloatingPointBitpattern(T) 1707 if (isFloatingPoint!T) 1708 { 1709 static if (T.mant_dig <= 64) 1710 { 1711 ulong mantissa; 1712 } 1713 else 1714 { 1715 ulong mantissa_lsb; 1716 ulong mantissa_msb; 1717 } 1718 1719 int exponent; 1720 bool negative; 1721 } 1722 1723 FloatingPointBitpattern!T extractBitpattern(T)(const(T) value) @trusted 1724 if (isFloatingPoint!T) 1725 { 1726 import std.math : floatTraits, RealFormat; 1727 1728 T val = value; 1729 FloatingPointBitpattern!T ret; 1730 1731 alias F = floatTraits!T; 1732 static if (F.realFormat == RealFormat.ieeeExtended) 1733 { 1734 if (__ctfe) 1735 { 1736 import core.math : fabs, ldexp; 1737 import std.math.rounding : floor; 1738 import std.math.traits : isInfinity, isNaN, signbit; 1739 import std.math.exponential : log2; 1740 1741 if (isNaN(val) || isInfinity(val)) 1742 ret.exponent = 32767; 1743 else if (fabs(val) < real.min_normal) 1744 ret.exponent = 0; 1745 else if (fabs(val) >= nextUp(real.max / 2)) 1746 ret.exponent = 32766; 1747 else 1748 ret.exponent = cast(int) (val.fabs.log2.floor() + 16383); 1749 1750 if (ret.exponent == 32767) 1751 { 1752 // NaN or infinity 1753 ret.mantissa = isNaN(val) ? ((1L << 63) - 1) : 0; 1754 } 1755 else 1756 { 1757 auto delta = 16382 + 64 // bias + bits of ulong 1758 - (ret.exponent == 0 ? 1 : ret.exponent); // -1 in case of subnormals 1759 val = ldexp(val, delta); // val *= 2^^delta 1760 1761 ulong tmp = cast(ulong) fabs(val); 1762 if (ret.exponent != 32767 && ret.exponent > 0 && tmp <= ulong.max / 2) 1763 { 1764 // correction, due to log2(val) being rounded up: 1765 ret.exponent--; 1766 val *= 2; 1767 tmp = cast(ulong) fabs(val); 1768 } 1769 1770 ret.mantissa = tmp & long.max; 1771 } 1772 1773 ret.negative = (signbit(val) == 1); 1774 } 1775 else 1776 { 1777 ushort* vs = cast(ushort*) &val; 1778 ret.mantissa = (cast(ulong*) vs)[0] & long.max; 1779 ret.exponent = vs[4] & short.max; 1780 ret.negative = (vs[4] >> 15) & 1; 1781 } 1782 } 1783 else 1784 { 1785 static if (F.realFormat == RealFormat.ieeeSingle) 1786 { 1787 ulong ival = *cast(uint*) &val; 1788 } 1789 else static if (F.realFormat == RealFormat.ieeeDouble) 1790 { 1791 ulong ival = *cast(ulong*) &val; 1792 } 1793 else 1794 { 1795 static assert(false, "Floating point type `" ~ F.realFormat ~ "` not supported."); 1796 } 1797 1798 import std.math.exponential : log2; 1799 enum log2_max_exp = cast(int) log2(T(T.max_exp)); 1800 1801 ret.mantissa = ival & ((1L << (T.mant_dig - 1)) - 1); 1802 ret.exponent = (ival >> (T.mant_dig - 1)) & ((1L << (log2_max_exp + 1)) - 1); 1803 ret.negative = (ival >> (T.mant_dig + log2_max_exp)) & 1; 1804 } 1805 1806 // add leading 1 for normalized values and correct exponent for denormalied values 1807 if (ret.exponent != 0 && ret.exponent != 2 * T.max_exp - 1) 1808 ret.mantissa |= 1L << (T.mant_dig - 1); 1809 else if (ret.exponent == 0) 1810 ret.exponent = 1; 1811 1812 ret.exponent -= T.max_exp - 1; 1813 1814 return ret; 1815 } 1816 1817 @safe pure unittest 1818 { 1819 float f = 1.0f; 1820 auto bp = extractBitpattern(f); 1821 assert(bp.mantissa == 0x80_0000); 1822 assert(bp.exponent == 0); 1823 assert(bp.negative == false); 1824 1825 f = float.max; 1826 bp = extractBitpattern(f); 1827 assert(bp.mantissa == 0xff_ffff); 1828 assert(bp.exponent == 127); 1829 assert(bp.negative == false); 1830 1831 f = -1.5432e-17f; 1832 bp = extractBitpattern(f); 1833 assert(bp.mantissa == 0x8e_55c8); 1834 assert(bp.exponent == -56); 1835 assert(bp.negative == true); 1836 1837 // using double literal due to https://issues.dlang.org/show_bug.cgi?id=20361 1838 f = 2.3822073893521890206e-44; 1839 bp = extractBitpattern(f); 1840 assert(bp.mantissa == 0x00_0011); 1841 assert(bp.exponent == -126); 1842 assert(bp.negative == false); 1843 1844 f = -float.infinity; 1845 bp = extractBitpattern(f); 1846 assert(bp.mantissa == 0); 1847 assert(bp.exponent == 128); 1848 assert(bp.negative == true); 1849 1850 f = float.nan; 1851 bp = extractBitpattern(f); 1852 assert(bp.mantissa != 0); // we don't guarantee payloads 1853 assert(bp.exponent == 128); 1854 assert(bp.negative == false); 1855 } 1856 1857 @safe pure unittest 1858 { 1859 double d = 1.0; 1860 auto bp = extractBitpattern(d); 1861 assert(bp.mantissa == 0x10_0000_0000_0000L); 1862 assert(bp.exponent == 0); 1863 assert(bp.negative == false); 1864 1865 d = double.max; 1866 bp = extractBitpattern(d); 1867 assert(bp.mantissa == 0x1f_ffff_ffff_ffffL); 1868 assert(bp.exponent == 1023); 1869 assert(bp.negative == false); 1870 1871 d = -1.5432e-222; 1872 bp = extractBitpattern(d); 1873 assert(bp.mantissa == 0x11_d9b6_a401_3b04L); 1874 assert(bp.exponent == -737); 1875 assert(bp.negative == true); 1876 1877 d = 0.0.nextUp; 1878 bp = extractBitpattern(d); 1879 assert(bp.mantissa == 0x00_0000_0000_0001L); 1880 assert(bp.exponent == -1022); 1881 assert(bp.negative == false); 1882 1883 d = -double.infinity; 1884 bp = extractBitpattern(d); 1885 assert(bp.mantissa == 0); 1886 assert(bp.exponent == 1024); 1887 assert(bp.negative == true); 1888 1889 d = double.nan; 1890 bp = extractBitpattern(d); 1891 assert(bp.mantissa != 0); // we don't guarantee payloads 1892 assert(bp.exponent == 1024); 1893 assert(bp.negative == false); 1894 } 1895 1896 @safe pure unittest 1897 { 1898 import std.math : floatTraits, RealFormat; 1899 1900 alias F = floatTraits!real; 1901 static if (F.realFormat == RealFormat.ieeeExtended) 1902 { 1903 real r = 1.0L; 1904 auto bp = extractBitpattern(r); 1905 assert(bp.mantissa == 0x8000_0000_0000_0000L); 1906 assert(bp.exponent == 0); 1907 assert(bp.negative == false); 1908 1909 r = real.max; 1910 bp = extractBitpattern(r); 1911 assert(bp.mantissa == 0xffff_ffff_ffff_ffffL); 1912 assert(bp.exponent == 16383); 1913 assert(bp.negative == false); 1914 1915 r = -1.5432e-3333L; 1916 bp = extractBitpattern(r); 1917 assert(bp.mantissa == 0xc768_a2c7_a616_cc22L); 1918 assert(bp.exponent == -11072); 1919 assert(bp.negative == true); 1920 1921 r = 0.0L.nextUp; 1922 bp = extractBitpattern(r); 1923 assert(bp.mantissa == 0x0000_0000_0000_0001L); 1924 assert(bp.exponent == -16382); 1925 assert(bp.negative == false); 1926 1927 r = -float.infinity; 1928 bp = extractBitpattern(r); 1929 assert(bp.mantissa == 0); 1930 assert(bp.exponent == 16384); 1931 assert(bp.negative == true); 1932 1933 r = float.nan; 1934 bp = extractBitpattern(r); 1935 assert(bp.mantissa != 0); // we don't guarantee payloads 1936 assert(bp.exponent == 16384); 1937 assert(bp.negative == false); 1938 1939 r = nextDown(0x1p+16383L); 1940 bp = extractBitpattern(r); 1941 assert(bp.mantissa == 0xffff_ffff_ffff_ffffL); 1942 assert(bp.exponent == 16382); 1943 assert(bp.negative == false); 1944 } 1945 } 1946 1947 @safe pure unittest 1948 { 1949 import std.math : floatTraits, RealFormat; 1950 import std.math.exponential : log2; 1951 1952 alias F = floatTraits!real; 1953 1954 // log2 is broken for x87-reals on some computers in CTFE 1955 // the following test excludes these computers from the test 1956 // (https://issues.dlang.org/show_bug.cgi?id=21757) 1957 enum test = cast(int) log2(3.05e2312L); 1958 static if (F.realFormat == RealFormat.ieeeExtended && test == 7681) 1959 { 1960 enum r1 = 1.0L; 1961 enum bp1 = extractBitpattern(r1); 1962 static assert(bp1.mantissa == 0x8000_0000_0000_0000L); 1963 static assert(bp1.exponent == 0); 1964 static assert(bp1.negative == false); 1965 1966 enum r2 = real.max; 1967 enum bp2 = extractBitpattern(r2); 1968 static assert(bp2.mantissa == 0xffff_ffff_ffff_ffffL); 1969 static assert(bp2.exponent == 16383); 1970 static assert(bp2.negative == false); 1971 1972 enum r3 = -1.5432e-3333L; 1973 enum bp3 = extractBitpattern(r3); 1974 static assert(bp3.mantissa == 0xc768_a2c7_a616_cc22L); 1975 static assert(bp3.exponent == -11072); 1976 static assert(bp3.negative == true); 1977 1978 enum r4 = 0.0L.nextUp; 1979 enum bp4 = extractBitpattern(r4); 1980 static assert(bp4.mantissa == 0x0000_0000_0000_0001L); 1981 static assert(bp4.exponent == -16382); 1982 static assert(bp4.negative == false); 1983 1984 enum r5 = -real.infinity; 1985 enum bp5 = extractBitpattern(r5); 1986 static assert(bp5.mantissa == 0); 1987 static assert(bp5.exponent == 16384); 1988 static assert(bp5.negative == true); 1989 1990 enum r6 = real.nan; 1991 enum bp6 = extractBitpattern(r6); 1992 static assert(bp6.mantissa != 0); // we don't guarantee payloads 1993 static assert(bp6.exponent == 16384); 1994 static assert(bp6.negative == false); 1995 1996 enum r7 = nextDown(0x1p+16383L); 1997 enum bp7 = extractBitpattern(r7); 1998 static assert(bp7.mantissa == 0xffff_ffff_ffff_ffffL); 1999 static assert(bp7.exponent == 16382); 2000 static assert(bp7.negative == false); 2001 } 2002 }