1// Written in the D programming language. 2 3/** 4This is a submodule of $(MREF std, math). 5 6It contains several functions for work with floating point numbers. 7 8Copyright: Copyright The D Language Foundation 2000 - 2011. 9License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). 10Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston, 11 Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger 12Source: $(PHOBOSSRC std/math/operations.d) 13 14Macros: 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 27module std.math.operations; 28 29import 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 */ 45real 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 */ 160ulong 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 242debug(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 */ 284real 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 */ 422double 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 */ 449float 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 */ 500real nextDown(real x) @safe pure nothrow @nogc 501{ 502 return -nextUp(-x); 503} 504 505/** ditto */ 506double nextDown(double x) @safe pure nothrow @nogc 507{ 508 return -nextUp(-x); 509} 510 511/** ditto */ 512float 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 */ 664T 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 */ 749real 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 */ 774F fmax(F)(const F x, const F y) @safe pure nothrow @nogc 775if (__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 */ 806F fmin(F)(const F x, const F y) @safe pure nothrow @nogc 807if (__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 */ 837pragma(inline, true) 838real 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 */ 865int feqrel(X)(const X x, const X y) @trusted pure nothrow @nogc 866if (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 */ 1039deprecated("approxEqual will be removed in 2.106.0. Please use isClose instead.") 1040bool 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 1109deprecated @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 1124deprecated @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 1137deprecated @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 1147deprecated @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 1157deprecated @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 1172deprecated @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 1186deprecated @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 */ 1225bool 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)) && is(typeof(rhs.infinity))) 1291 { 1292 if (lhs == lhs.infinity || rhs == rhs.infinity || 1293 lhs == -lhs.infinity || rhs == -rhs.infinity) return false; 1294 } 1295 1296 import std.math.algebraic : abs; 1297 1298 auto diff = abs(lhs - rhs); 1299 1300 return diff <= maxRelDiff*abs(lhs) 1301 || diff <= maxRelDiff*abs(rhs) 1302 || diff <= maxAbsDiff; 1303 } 1304} 1305 1306/// 1307@safe pure nothrow @nogc unittest 1308{ 1309 assert(isClose(1.0,0.999_999_999)); 1310 assert(isClose(0.001, 0.000_999_999_999)); 1311 assert(isClose(1_000_000_000.0,999_999_999.0)); 1312 1313 assert(isClose(17.123_456_789, 17.123_456_78)); 1314 assert(!isClose(17.123_456_789, 17.123_45)); 1315 1316 // use explicit 3rd parameter for less (or more) accuracy 1317 assert(isClose(17.123_456_789, 17.123_45, 1e-6)); 1318 assert(!isClose(17.123_456_789, 17.123_45, 1e-7)); 1319 1320 // use 4th parameter when comparing close to zero 1321 assert(!isClose(1e-100, 0.0)); 1322 assert(isClose(1e-100, 0.0, 0.0, 1e-90)); 1323 assert(!isClose(1e-10, -1e-10)); 1324 assert(isClose(1e-10, -1e-10, 0.0, 1e-9)); 1325 assert(!isClose(1e-300, 1e-298)); 1326 assert(isClose(1e-300, 1e-298, 0.0, 1e-200)); 1327 1328 // different default limits for different floating point types 1329 assert(isClose(1.0f, 0.999_99f)); 1330 assert(!isClose(1.0, 0.999_99)); 1331 static if (real.sizeof > double.sizeof) 1332 assert(!isClose(1.0L, 0.999_999_999L)); 1333} 1334 1335/// 1336@safe pure nothrow unittest 1337{ 1338 assert(isClose([1.0, 2.0, 3.0], [0.999_999_999, 2.000_000_001, 3.0])); 1339 assert(!isClose([1.0, 2.0], [0.999_999_999, 2.000_000_001, 3.0])); 1340 assert(!isClose([1.0, 2.0, 3.0], [0.999_999_999, 2.000_000_001])); 1341 1342 assert(isClose([2.0, 1.999_999_999, 2.000_000_001], 2.0)); 1343 assert(isClose(2.0, [2.0, 1.999_999_999, 2.000_000_001])); 1344} 1345 1346@safe pure nothrow unittest 1347{ 1348 assert(!isClose([1.0, 2.0, 3.0], [0.999_999_999, 3.0, 3.0])); 1349 assert(!isClose([2.0, 1.999_999, 2.000_000_001], 2.0)); 1350 assert(!isClose(2.0, [2.0, 1.999_999_999, 2.000_000_999])); 1351} 1352 1353@safe pure nothrow @nogc unittest 1354{ 1355 immutable a = 1.00001f; 1356 const b = 1.000019; 1357 assert(isClose(a,b)); 1358 1359 assert(isClose(1.00001f,1.000019f)); 1360 assert(isClose(1.00001f,1.000019)); 1361 assert(isClose(1.00001,1.000019f)); 1362 assert(!isClose(1.00001,1.000019)); 1363 1364 real a1 = 1e-300L; 1365 real a2 = a1.nextUp; 1366 assert(isClose(a1,a2)); 1367} 1368 1369@safe pure nothrow unittest 1370{ 1371 float[] arr1 = [ 1.0, 2.0, 3.0 ]; 1372 double[] arr2 = [ 1.00001, 1.99999, 3 ]; 1373 assert(isClose(arr1, arr2)); 1374} 1375 1376@safe pure nothrow @nogc unittest 1377{ 1378 assert(!isClose(1000.0,1010.0)); 1379 assert(!isClose(9_090_000_000.0,9_000_000_000.0)); 1380 assert(isClose(0.0,1e30,1.0)); 1381 assert(!isClose(0.00001,1e-30)); 1382 assert(!isClose(-1e-30,1e-30,1e-2,0.0)); 1383} 1384 1385@safe pure nothrow @nogc unittest 1386{ 1387 assert(!isClose(3, 0)); 1388 assert(isClose(3, 3)); 1389 assert(isClose(3.0, 3)); 1390 assert(isClose(3, 3.0)); 1391 1392 assert(isClose(0.0,0.0)); 1393 assert(isClose(-0.0,0.0)); 1394 assert(isClose(0.0f,0.0)); 1395} 1396 1397@safe pure nothrow @nogc unittest 1398{ 1399 real num = real.infinity; 1400 assert(num == real.infinity); 1401 assert(isClose(num, real.infinity)); 1402 num = -real.infinity; 1403 assert(num == -real.infinity); 1404 assert(isClose(num, -real.infinity)); 1405 1406 assert(!isClose(1,real.nan)); 1407 assert(!isClose(real.nan,real.max)); 1408 assert(!isClose(real.nan,real.nan)); 1409} 1410 1411@safe pure nothrow @nogc unittest 1412{ 1413 assert(isClose!(real[],real[],real)([],[])); 1414 assert(isClose(cast(real[])[],cast(real[])[])); 1415} 1416 1417@safe pure nothrow @nogc unittest 1418{ 1419 import std.conv : to; 1420 1421 float f = 31.79f; 1422 double d = 31.79; 1423 double f2d = f.to!double; 1424 1425 assert(isClose(f,f2d)); 1426 assert(!isClose(d,f2d)); 1427} 1428 1429@safe pure nothrow @nogc unittest 1430{ 1431 import std.conv : to; 1432 1433 double d = 31.79; 1434 float f = d.to!float; 1435 double f2d = f.to!double; 1436 1437 assert(isClose(f,f2d)); 1438 assert(!isClose(d,f2d)); 1439 assert(isClose(d,f2d,1e-4)); 1440} 1441 1442package(std.math) template CommonDefaultFor(T,U) 1443{ 1444 import std.algorithm.comparison : min; 1445 1446 alias baseT = FloatingPointBaseType!T; 1447 alias baseU = FloatingPointBaseType!U; 1448 1449 enum CommonType!(baseT, baseU) CommonDefaultFor = 10.0L ^^ -((min(baseT.dig, baseU.dig) + 1) / 2 + 1); 1450} 1451 1452private template FloatingPointBaseType(T) 1453{ 1454 import std.range.primitives : ElementType; 1455 static if (isFloatingPoint!T) 1456 { 1457 alias FloatingPointBaseType = Unqual!T; 1458 } 1459 else static if (isFloatingPoint!(ElementType!(Unqual!T))) 1460 { 1461 alias FloatingPointBaseType = Unqual!(ElementType!(Unqual!T)); 1462 } 1463 else 1464 { 1465 alias FloatingPointBaseType = real; 1466 } 1467} 1468 1469/*********************************** 1470 * Defines a total order on all floating-point numbers. 1471 * 1472 * The order is defined as follows: 1473 * $(UL 1474 * $(LI All numbers in [-$(INFIN), +$(INFIN)] are ordered 1475 * the same way as by built-in comparison, with the exception of 1476 * -0.0, which is less than +0.0;) 1477 * $(LI If the sign bit is set (that is, it's 'negative'), $(NAN) is less 1478 * than any number; if the sign bit is not set (it is 'positive'), 1479 * $(NAN) is greater than any number;) 1480 * $(LI $(NAN)s of the same sign are ordered by the payload ('negative' 1481 * ones - in reverse order).) 1482 * ) 1483 * 1484 * Returns: 1485 * negative value if `x` precedes `y` in the order specified above; 1486 * 0 if `x` and `y` are identical, and positive value otherwise. 1487 * 1488 * See_Also: 1489 * $(MYREF isIdentical) 1490 * Standards: Conforms to IEEE 754-2008 1491 */ 1492int cmp(T)(const(T) x, const(T) y) @nogc @trusted pure nothrow 1493if (isFloatingPoint!T) 1494{ 1495 import std.math : floatTraits, RealFormat; 1496 1497 alias F = floatTraits!T; 1498 1499 static if (F.realFormat == RealFormat.ieeeSingle 1500 || F.realFormat == RealFormat.ieeeDouble) 1501 { 1502 static if (T.sizeof == 4) 1503 alias UInt = uint; 1504 else 1505 alias UInt = ulong; 1506 1507 union Repainter 1508 { 1509 T number; 1510 UInt bits; 1511 } 1512 1513 enum msb = ~(UInt.max >>> 1); 1514 1515 import std.typecons : Tuple; 1516 Tuple!(Repainter, Repainter) vars = void; 1517 vars[0].number = x; 1518 vars[1].number = y; 1519 1520 foreach (ref var; vars) 1521 if (var.bits & msb) 1522 var.bits = ~var.bits; 1523 else 1524 var.bits |= msb; 1525 1526 if (vars[0].bits < vars[1].bits) 1527 return -1; 1528 else if (vars[0].bits > vars[1].bits) 1529 return 1; 1530 else 1531 return 0; 1532 } 1533 else static if (F.realFormat == RealFormat.ieeeExtended53 1534 || F.realFormat == RealFormat.ieeeExtended 1535 || F.realFormat == RealFormat.ieeeQuadruple) 1536 { 1537 static if (F.realFormat == RealFormat.ieeeQuadruple) 1538 alias RemT = ulong; 1539 else 1540 alias RemT = ushort; 1541 1542 struct Bits 1543 { 1544 ulong bulk; 1545 RemT rem; 1546 } 1547 1548 union Repainter 1549 { 1550 T number; 1551 Bits bits; 1552 ubyte[T.sizeof] bytes; 1553 } 1554 1555 import std.typecons : Tuple; 1556 Tuple!(Repainter, Repainter) vars = void; 1557 vars[0].number = x; 1558 vars[1].number = y; 1559 1560 foreach (ref var; vars) 1561 if (var.bytes[F.SIGNPOS_BYTE] & 0x80) 1562 { 1563 var.bits.bulk = ~var.bits.bulk; 1564 var.bits.rem = cast(typeof(var.bits.rem))(-1 - var.bits.rem); // ~var.bits.rem 1565 } 1566 else 1567 { 1568 var.bytes[F.SIGNPOS_BYTE] |= 0x80; 1569 } 1570 1571 version (LittleEndian) 1572 { 1573 if (vars[0].bits.rem < vars[1].bits.rem) 1574 return -1; 1575 else if (vars[0].bits.rem > vars[1].bits.rem) 1576 return 1; 1577 else if (vars[0].bits.bulk < vars[1].bits.bulk) 1578 return -1; 1579 else if (vars[0].bits.bulk > vars[1].bits.bulk) 1580 return 1; 1581 else 1582 return 0; 1583 } 1584 else 1585 { 1586 if (vars[0].bits.bulk < vars[1].bits.bulk) 1587 return -1; 1588 else if (vars[0].bits.bulk > vars[1].bits.bulk) 1589 return 1; 1590 else if (vars[0].bits.rem < vars[1].bits.rem) 1591 return -1; 1592 else if (vars[0].bits.rem > vars[1].bits.rem) 1593 return 1; 1594 else 1595 return 0; 1596 } 1597 } 1598 else 1599 { 1600 // IBM Extended doubledouble does not follow the general 1601 // sign-exponent-significand layout, so has to be handled generically 1602 1603 import std.math.traits : signbit, isNaN; 1604 1605 const int xSign = signbit(x), 1606 ySign = signbit(y); 1607 1608 if (xSign == 1 && ySign == 1) 1609 return cmp(-y, -x); 1610 else if (xSign == 1) 1611 return -1; 1612 else if (ySign == 1) 1613 return 1; 1614 else if (x < y) 1615 return -1; 1616 else if (x == y) 1617 return 0; 1618 else if (x > y) 1619 return 1; 1620 else if (isNaN(x) && !isNaN(y)) 1621 return 1; 1622 else if (isNaN(y) && !isNaN(x)) 1623 return -1; 1624 else if (getNaNPayload(x) < getNaNPayload(y)) 1625 return -1; 1626 else if (getNaNPayload(x) > getNaNPayload(y)) 1627 return 1; 1628 else 1629 return 0; 1630 } 1631} 1632 1633/// Most numbers are ordered naturally. 1634@safe unittest 1635{ 1636 assert(cmp(-double.infinity, -double.max) < 0); 1637 assert(cmp(-double.max, -100.0) < 0); 1638 assert(cmp(-100.0, -0.5) < 0); 1639 assert(cmp(-0.5, 0.0) < 0); 1640 assert(cmp(0.0, 0.5) < 0); 1641 assert(cmp(0.5, 100.0) < 0); 1642 assert(cmp(100.0, double.max) < 0); 1643 assert(cmp(double.max, double.infinity) < 0); 1644 1645 assert(cmp(1.0, 1.0) == 0); 1646} 1647 1648/// Positive and negative zeroes are distinct. 1649@safe unittest 1650{ 1651 assert(cmp(-0.0, +0.0) < 0); 1652 assert(cmp(+0.0, -0.0) > 0); 1653} 1654 1655/// Depending on the sign, $(NAN)s go to either end of the spectrum. 1656@safe unittest 1657{ 1658 assert(cmp(-double.nan, -double.infinity) < 0); 1659 assert(cmp(double.infinity, double.nan) < 0); 1660 assert(cmp(-double.nan, double.nan) < 0); 1661} 1662 1663/// $(NAN)s of the same sign are ordered by the payload. 1664@safe unittest 1665{ 1666 assert(cmp(NaN(10), NaN(20)) < 0); 1667 assert(cmp(-NaN(20), -NaN(10)) < 0); 1668} 1669 1670@safe unittest 1671{ 1672 import std.meta : AliasSeq; 1673 static foreach (T; AliasSeq!(float, double, real)) 1674 {{ 1675 T[] values = [-cast(T) NaN(20), -cast(T) NaN(10), -T.nan, -T.infinity, 1676 -T.max, -T.max / 2, T(-16.0), T(-1.0).nextDown, 1677 T(-1.0), T(-1.0).nextUp, 1678 T(-0.5), -T.min_normal, (-T.min_normal).nextUp, 1679 -2 * T.min_normal * T.epsilon, 1680 -T.min_normal * T.epsilon, 1681 T(-0.0), T(0.0), 1682 T.min_normal * T.epsilon, 1683 2 * T.min_normal * T.epsilon, 1684 T.min_normal.nextDown, T.min_normal, T(0.5), 1685 T(1.0).nextDown, T(1.0), 1686 T(1.0).nextUp, T(16.0), T.max / 2, T.max, 1687 T.infinity, T.nan, cast(T) NaN(10), cast(T) NaN(20)]; 1688 1689 foreach (i, x; values) 1690 { 1691 foreach (y; values[i + 1 .. $]) 1692 { 1693 assert(cmp(x, y) < 0); 1694 assert(cmp(y, x) > 0); 1695 } 1696 assert(cmp(x, x) == 0); 1697 } 1698 }} 1699} 1700 1701package(std): // not yet public 1702 1703struct FloatingPointBitpattern(T) 1704if (isFloatingPoint!T) 1705{ 1706 static if (T.mant_dig <= 64) 1707 { 1708 ulong mantissa; 1709 } 1710 else 1711 { 1712 ulong mantissa_lsb; 1713 ulong mantissa_msb; 1714 } 1715 1716 int exponent; 1717 bool negative; 1718} 1719 1720FloatingPointBitpattern!T extractBitpattern(T)(const(T) value) @trusted 1721if (isFloatingPoint!T) 1722{ 1723 import std.math : floatTraits, RealFormat; 1724 1725 T val = value; 1726 FloatingPointBitpattern!T ret; 1727 1728 alias F = floatTraits!T; 1729 static if (F.realFormat == RealFormat.ieeeExtended) 1730 { 1731 if (__ctfe) 1732 { 1733 import core.math : fabs, ldexp; 1734 import std.math.rounding : floor; 1735 import std.math.traits : isInfinity, isNaN, signbit; 1736 import std.math.exponential : log2; 1737 1738 if (isNaN(val) || isInfinity(val)) 1739 ret.exponent = 32767; 1740 else if (fabs(val) < real.min_normal) 1741 ret.exponent = 0; 1742 else if (fabs(val) >= nextUp(real.max / 2)) 1743 ret.exponent = 32766; 1744 else 1745 ret.exponent = cast(int) (val.fabs.log2.floor() + 16383); 1746 1747 if (ret.exponent == 32767) 1748 { 1749 // NaN or infinity 1750 ret.mantissa = isNaN(val) ? ((1L << 63) - 1) : 0; 1751 } 1752 else 1753 { 1754 auto delta = 16382 + 64 // bias + bits of ulong 1755 - (ret.exponent == 0 ? 1 : ret.exponent); // -1 in case of subnormals 1756 val = ldexp(val, delta); // val *= 2^^delta 1757 1758 ulong tmp = cast(ulong) fabs(val); 1759 if (ret.exponent != 32767 && ret.exponent > 0 && tmp <= ulong.max / 2) 1760 { 1761 // correction, due to log2(val) being rounded up: 1762 ret.exponent--; 1763 val *= 2; 1764 tmp = cast(ulong) fabs(val); 1765 } 1766 1767 ret.mantissa = tmp & long.max; 1768 } 1769 1770 ret.negative = (signbit(val) == 1); 1771 } 1772 else 1773 { 1774 ushort* vs = cast(ushort*) &val; 1775 ret.mantissa = (cast(ulong*) vs)[0] & long.max; 1776 ret.exponent = vs[4] & short.max; 1777 ret.negative = (vs[4] >> 15) & 1; 1778 } 1779 } 1780 else 1781 { 1782 static if (F.realFormat == RealFormat.ieeeSingle) 1783 { 1784 ulong ival = *cast(uint*) &val; 1785 } 1786 else static if (F.realFormat == RealFormat.ieeeDouble) 1787 { 1788 ulong ival = *cast(ulong*) &val; 1789 } 1790 else 1791 { 1792 static assert(false, "Floating point type `" ~ F.realFormat ~ "` not supported."); 1793 } 1794 1795 import std.math.exponential : log2; 1796 enum log2_max_exp = cast(int) log2(T.max_exp); 1797 1798 ret.mantissa = ival & ((1L << (T.mant_dig - 1)) - 1); 1799 ret.exponent = (ival >> (T.mant_dig - 1)) & ((1L << (log2_max_exp + 1)) - 1); 1800 ret.negative = (ival >> (T.mant_dig + log2_max_exp)) & 1; 1801 } 1802 1803 // add leading 1 for normalized values and correct exponent for denormalied values 1804 if (ret.exponent != 0 && ret.exponent != 2 * T.max_exp - 1) 1805 ret.mantissa |= 1L << (T.mant_dig - 1); 1806 else if (ret.exponent == 0) 1807 ret.exponent = 1; 1808 1809 ret.exponent -= T.max_exp - 1; 1810 1811 return ret; 1812} 1813 1814@safe pure unittest 1815{ 1816 float f = 1.0f; 1817 auto bp = extractBitpattern(f); 1818 assert(bp.mantissa == 0x80_0000); 1819 assert(bp.exponent == 0); 1820 assert(bp.negative == false); 1821 1822 f = float.max; 1823 bp = extractBitpattern(f); 1824 assert(bp.mantissa == 0xff_ffff); 1825 assert(bp.exponent == 127); 1826 assert(bp.negative == false); 1827 1828 f = -1.5432e-17f; 1829 bp = extractBitpattern(f); 1830 assert(bp.mantissa == 0x8e_55c8); 1831 assert(bp.exponent == -56); 1832 assert(bp.negative == true); 1833 1834 // using double literal due to https://issues.dlang.org/show_bug.cgi?id=20361 1835 f = 2.3822073893521890206e-44; 1836 bp = extractBitpattern(f); 1837 assert(bp.mantissa == 0x00_0011); 1838 assert(bp.exponent == -126); 1839 assert(bp.negative == false); 1840 1841 f = -float.infinity; 1842 bp = extractBitpattern(f); 1843 assert(bp.mantissa == 0); 1844 assert(bp.exponent == 128); 1845 assert(bp.negative == true); 1846 1847 f = float.nan; 1848 bp = extractBitpattern(f); 1849 assert(bp.mantissa != 0); // we don't guarantee payloads 1850 assert(bp.exponent == 128); 1851 assert(bp.negative == false); 1852} 1853 1854@safe pure unittest 1855{ 1856 double d = 1.0; 1857 auto bp = extractBitpattern(d); 1858 assert(bp.mantissa == 0x10_0000_0000_0000L); 1859 assert(bp.exponent == 0); 1860 assert(bp.negative == false); 1861 1862 d = double.max; 1863 bp = extractBitpattern(d); 1864 assert(bp.mantissa == 0x1f_ffff_ffff_ffffL); 1865 assert(bp.exponent == 1023); 1866 assert(bp.negative == false); 1867 1868 d = -1.5432e-222; 1869 bp = extractBitpattern(d); 1870 assert(bp.mantissa == 0x11_d9b6_a401_3b04L); 1871 assert(bp.exponent == -737); 1872 assert(bp.negative == true); 1873 1874 d = 0.0.nextUp; 1875 bp = extractBitpattern(d); 1876 assert(bp.mantissa == 0x00_0000_0000_0001L); 1877 assert(bp.exponent == -1022); 1878 assert(bp.negative == false); 1879 1880 d = -double.infinity; 1881 bp = extractBitpattern(d); 1882 assert(bp.mantissa == 0); 1883 assert(bp.exponent == 1024); 1884 assert(bp.negative == true); 1885 1886 d = double.nan; 1887 bp = extractBitpattern(d); 1888 assert(bp.mantissa != 0); // we don't guarantee payloads 1889 assert(bp.exponent == 1024); 1890 assert(bp.negative == false); 1891} 1892 1893@safe pure unittest 1894{ 1895 import std.math : floatTraits, RealFormat; 1896 1897 alias F = floatTraits!real; 1898 static if (F.realFormat == RealFormat.ieeeExtended) 1899 { 1900 real r = 1.0L; 1901 auto bp = extractBitpattern(r); 1902 assert(bp.mantissa == 0x8000_0000_0000_0000L); 1903 assert(bp.exponent == 0); 1904 assert(bp.negative == false); 1905 1906 r = real.max; 1907 bp = extractBitpattern(r); 1908 assert(bp.mantissa == 0xffff_ffff_ffff_ffffL); 1909 assert(bp.exponent == 16383); 1910 assert(bp.negative == false); 1911 1912 r = -1.5432e-3333L; 1913 bp = extractBitpattern(r); 1914 assert(bp.mantissa == 0xc768_a2c7_a616_cc22L); 1915 assert(bp.exponent == -11072); 1916 assert(bp.negative == true); 1917 1918 r = 0.0L.nextUp; 1919 bp = extractBitpattern(r); 1920 assert(bp.mantissa == 0x0000_0000_0000_0001L); 1921 assert(bp.exponent == -16382); 1922 assert(bp.negative == false); 1923 1924 r = -float.infinity; 1925 bp = extractBitpattern(r); 1926 assert(bp.mantissa == 0); 1927 assert(bp.exponent == 16384); 1928 assert(bp.negative == true); 1929 1930 r = float.nan; 1931 bp = extractBitpattern(r); 1932 assert(bp.mantissa != 0); // we don't guarantee payloads 1933 assert(bp.exponent == 16384); 1934 assert(bp.negative == false); 1935 1936 r = nextDown(0x1p+16383L); 1937 bp = extractBitpattern(r); 1938 assert(bp.mantissa == 0xffff_ffff_ffff_ffffL); 1939 assert(bp.exponent == 16382); 1940 assert(bp.negative == false); 1941 } 1942} 1943 1944@safe pure unittest 1945{ 1946 import std.math : floatTraits, RealFormat; 1947 import std.math.exponential : log2; 1948 1949 alias F = floatTraits!real; 1950 1951 // log2 is broken for x87-reals on some computers in CTFE 1952 // the following test excludes these computers from the test 1953 // (issue 21757) 1954 enum test = cast(int) log2(3.05e2312L); 1955 static if (F.realFormat == RealFormat.ieeeExtended && test == 7681) 1956 { 1957 enum r1 = 1.0L; 1958 enum bp1 = extractBitpattern(r1); 1959 static assert(bp1.mantissa == 0x8000_0000_0000_0000L); 1960 static assert(bp1.exponent == 0); 1961 static assert(bp1.negative == false); 1962 1963 enum r2 = real.max; 1964 enum bp2 = extractBitpattern(r2); 1965 static assert(bp2.mantissa == 0xffff_ffff_ffff_ffffL); 1966 static assert(bp2.exponent == 16383); 1967 static assert(bp2.negative == false); 1968 1969 enum r3 = -1.5432e-3333L; 1970 enum bp3 = extractBitpattern(r3); 1971 static assert(bp3.mantissa == 0xc768_a2c7_a616_cc22L); 1972 static assert(bp3.exponent == -11072); 1973 static assert(bp3.negative == true); 1974 1975 enum r4 = 0.0L.nextUp; 1976 enum bp4 = extractBitpattern(r4); 1977 static assert(bp4.mantissa == 0x0000_0000_0000_0001L); 1978 static assert(bp4.exponent == -16382); 1979 static assert(bp4.negative == false); 1980 1981 enum r5 = -real.infinity; 1982 enum bp5 = extractBitpattern(r5); 1983 static assert(bp5.mantissa == 0); 1984 static assert(bp5.exponent == 16384); 1985 static assert(bp5.negative == true); 1986 1987 enum r6 = real.nan; 1988 enum bp6 = extractBitpattern(r6); 1989 static assert(bp6.mantissa != 0); // we don't guarantee payloads 1990 static assert(bp6.exponent == 16384); 1991 static assert(bp6.negative == false); 1992 1993 enum r7 = nextDown(0x1p+16383L); 1994 enum bp7 = extractBitpattern(r7); 1995 static assert(bp7.mantissa == 0xffff_ffff_ffff_ffffL); 1996 static assert(bp7.exponent == 16382); 1997 static assert(bp7.negative == false); 1998 } 1999} 2000