1 /** 2 * Low-level Mathematical Functions which take advantage of the IEEE754 ABI. 3 * 4 * Copyright: Portions Copyright (C) 2001-2005 Digital Mars. 5 * License: BSD style: $(LICENSE), Digital Mars. 6 * Authors: Don Clugston, Walter Bright, Sean Kelly 7 */ 8 /* Portions of this code were taken from Phobos std.math, which has the following 9 * copyright notice: 10 * 11 * Author: 12 * Walter Bright 13 * Copyright: 14 * Copyright (c) 2001-2005 by Digital Mars, 15 * All Rights Reserved, 16 * www.digitalmars.com 17 * License: 18 * This software is provided 'as-is', without any express or implied 19 * warranty. In no event will the authors be held liable for any damages 20 * arising from the use of this software. 21 * 22 * Permission is granted to anyone to use this software for any purpose, 23 * including commercial applications, and to alter it and redistribute it 24 * freely, subject to the following restrictions: 25 * 26 * <ul> 27 * <li> The origin of this software must not be misrepresented; you must not 28 * claim that you wrote the original software. If you use this software 29 * in a product, an acknowledgment in the product documentation would be 30 * appreciated but is not required. 31 * </li> 32 * <li> Altered source versions must be plainly marked as such, and must not 33 * be misrepresented as being the original software. 34 * </li> 35 * <li> This notice may not be removed or altered from any source 36 * distribution. 37 * </li> 38 * </ul> 39 */ 40 /** 41 * Macros: 42 * 43 * TABLE_SV = <table border=1 cellpadding=4 cellspacing=0> 44 * <caption>Special Values</caption> 45 * $0</table> 46 * SVH = $(TR $(TH $1) $(TH $2)) 47 * SV = $(TR $(TD $1) $(TD $2)) 48 * SVH3 = $(TR $(TH $1) $(TH $2) $(TH $3)) 49 * SV3 = $(TR $(TD $1) $(TD $2) $(TD $3)) 50 * NAN = $(RED NAN) 51 * PLUSMN = ± 52 * INFIN = ∞ 53 * PLUSMNINF = ±∞ 54 * PI = π 55 * LT = < 56 * GT = > 57 * SQRT = &radix; 58 * HALF = ½ 59 */ 60 module tango.math.IEEE; 61 62 version(GNU){ 63 // GDC is a filthy liar. It can't actually do inline asm. 64 } else version(TangoNoAsm) { 65 66 } else version(D_InlineAsm_X86) { 67 version = Naked_D_InlineAsm_X86; 68 } 69 70 version (X86){ 71 version = X86_Any; 72 } 73 74 version (X86_64){ 75 version = X86_Any; 76 } 77 78 version (Naked_D_InlineAsm_X86) { 79 // Don't include this extra dependency unless we need to. 80 debug(UnitTest) { 81 static import tango.stdc.math; 82 } 83 } else { 84 // Needed for cos(), sin(), tan() on GNU. 85 static import tango.stdc.math; 86 } 87 88 89 version(Windows) { 90 version(DigitalMars) { 91 version = DMDWindows; 92 } 93 } 94 95 // Standard Tango NaN payloads. 96 // NOTE: These values may change in future Tango releases 97 // The lowest three bits indicate the cause of the NaN: 98 // 0 = error other than those listed below: 99 // 1 = domain error 100 // 2 = singularity 101 // 3 = range 102 // 4-7 = reserved. 103 enum TANGO_NAN { 104 // General errors 105 DOMAIN_ERROR = 0x0101, 106 SINGULARITY = 0x0102, 107 RANGE_ERROR = 0x0103, 108 // NaNs created by functions in the basic library 109 TAN_DOMAIN = 0x1001, 110 POW_DOMAIN = 0x1021, 111 GAMMA_DOMAIN = 0x1101, 112 GAMMA_POLE = 0x1102, 113 SGNGAMMA = 0x1112, 114 BETA_DOMAIN = 0x1131, 115 // NaNs from statistical functions 116 NORMALDISTRIBUTION_INV_DOMAIN = 0x2001, 117 STUDENTSDDISTRIBUTION_DOMAIN = 0x2011 118 } 119 120 private: 121 /* Most of the functions depend on the format of the largest IEEE floating-point type. 122 * These code will differ depending on whether 'real' is 64, 80, or 128 bits, 123 * and whether it is a big-endian or little-endian architecture. 124 * Only five 'real' ABIs are currently supported: 125 * 64 bit Big-endian 'double' (eg PowerPC) 126 * 128 bit Big-endian 'quadruple' (eg SPARC) 127 * 64 bit Little-endian 'double' (eg x86-SSE2) 128 * 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium). 129 * 128 bit Little-endian 'quadruple' (not implemented on any known processor!) 130 * 131 * There is also an unsupported ABI which does not follow IEEE; several of its functions 132 * will generate run-time errors if used. 133 * 128 bit Big-endian 'doubledouble' (used by GDC <= 0.23 for PowerPC) 134 */ 135 136 version(LittleEndian) { 137 static assert(real.mant_dig == 53 || real.mant_dig==64 || real.mant_dig == 113, 138 "Only 64-bit, 80-bit, and 128-bit reals are supported for LittleEndian CPUs"); 139 } else { 140 static assert(real.mant_dig == 53 || real.mant_dig==106 || real.mant_dig == 113, 141 "Only 64-bit and 128-bit reals are supported for BigEndian CPUs. double-double reals have partial support"); 142 } 143 144 // Constants used for extracting the components of the representation. 145 // They supplement the built-in floating point properties. 146 template floatTraits(T) { 147 // EXPMASK is a ushort mask to select the exponent portion (without sign) 148 // SIGNMASK is a ushort mask to select the sign bit. 149 // EXPPOS_SHORT is the index of the exponent when represented as a ushort array. 150 // SIGNPOS_BYTE is the index of the sign when represented as a ubyte array. 151 // RECIP_EPSILON is the value such that (smallest_denormal) * RECIP_EPSILON == T.min 152 enum T RECIP_EPSILON = (1/T.epsilon); 153 154 static if (T.mant_dig == 24) { // float 155 enum : ushort { 156 EXPMASK = 0x7F80, 157 SIGNMASK = 0x8000, 158 EXPBIAS = 0x3F00 159 } 160 enum uint EXPMASK_INT = 0x7F80_0000; 161 enum uint MANTISSAMASK_INT = 0x007F_FFFF; 162 version(LittleEndian) { 163 enum EXPPOS_SHORT = 1; 164 } else { 165 enum EXPPOS_SHORT = 0; 166 } 167 } else static if (T.mant_dig==53) { // double, or real==double 168 enum : ushort { 169 EXPMASK = 0x7FF0, 170 SIGNMASK = 0x8000, 171 EXPBIAS = 0x3FE0 172 } 173 enum uint EXPMASK_INT = 0x7FF0_0000; 174 enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only 175 version(LittleEndian) { 176 enum EXPPOS_SHORT = 3; 177 enum SIGNPOS_BYTE = 7; 178 } else { 179 enum EXPPOS_SHORT = 0; 180 enum SIGNPOS_BYTE = 0; 181 } 182 } else static if (T.mant_dig==64) { // real80 183 enum : ushort { 184 EXPMASK = 0x7FFF, 185 SIGNMASK = 0x8000, 186 EXPBIAS = 0x3FFE 187 } 188 // enum ulong QUIETNANMASK = 0xC000_0000_0000_0000; // Converts a signaling NaN to a quiet NaN. 189 version(LittleEndian) { 190 enum EXPPOS_SHORT = 4; 191 enum SIGNPOS_BYTE = 9; 192 } else { 193 enum EXPPOS_SHORT = 0; 194 enum SIGNPOS_BYTE = 0; 195 } 196 } else static if (real.mant_dig==113){ // quadruple 197 enum : ushort { 198 EXPMASK = 0x7FFF, 199 SIGNMASK = 0x8000, 200 EXPBIAS = 0x3FFE 201 } 202 version(LittleEndian) { 203 enum EXPPOS_SHORT = 7; 204 enum SIGNPOS_BYTE = 15; 205 } else { 206 enum EXPPOS_SHORT = 0; 207 enum SIGNPOS_BYTE = 0; 208 } 209 } else static if (real.mant_dig==106) { // doubledouble 210 enum : ushort { 211 EXPMASK = 0x7FF0, 212 SIGNMASK = 0x8000 213 // EXPBIAS = 0x3FE0 214 } 215 // the exponent byte is not unique 216 version(LittleEndian) { 217 enum EXPPOS_SHORT = 7; // 3 is also an exp short 218 enum SIGNPOS_BYTE = 15; 219 } else { 220 enum EXPPOS_SHORT = 0; // 4 is also an exp short 221 enum SIGNPOS_BYTE = 0; 222 } 223 } 224 } 225 226 // These apply to all floating-point types 227 version(LittleEndian) { 228 enum MANTISSA_LSB = 0; 229 enum MANTISSA_MSB = 1; 230 } else { 231 enum MANTISSA_LSB = 1; 232 enum MANTISSA_MSB = 0; 233 } 234 235 public: 236 237 /** IEEE exception status flags 238 239 These flags indicate that an exceptional floating-point condition has occured. 240 They indicate that a NaN or an infinity has been generated, that a result 241 is inexact, or that a signalling NaN has been encountered. 242 The return values of the properties should be treated as booleans, although 243 each is returned as an int, for speed. 244 245 Example: 246 ---- 247 real a=3.5; 248 // Set all the flags to zero 249 resetIeeeFlags(); 250 assert(!ieeeFlags.divByZero); 251 // Perform a division by zero. 252 a/=0.0L; 253 assert(a==real.infinity); 254 assert(ieeeFlags.divByZero); 255 // Create a NaN 256 a*=0.0L; 257 assert(ieeeFlags.invalid); 258 assert(isNaN(a)); 259 260 // Check that calling func() has no effect on the 261 // status flags. 262 IeeeFlags f = ieeeFlags; 263 func(); 264 assert(ieeeFlags == f); 265 266 ---- 267 */ 268 struct IeeeFlags 269 { 270 private: 271 // The x87 FPU status register is 16 bits. 272 // The Pentium SSE2 status register is 32 bits. 273 int m_flags; 274 version (X86_Any) { 275 // Applies to both x87 status word (16 bits) and SSE2 status word(32 bits). 276 enum : int { 277 INEXACT_MASK = 0x20, 278 UNDERFLOW_MASK = 0x10, 279 OVERFLOW_MASK = 0x08, 280 DIVBYZERO_MASK = 0x04, 281 INVALID_MASK = 0x01 282 } 283 // Don't bother about denormals, they are not supported on most CPUs. 284 // DENORMAL_MASK = 0x02; 285 } else version (PPC) { 286 // PowerPC FPSCR is a 32-bit register. 287 enum : int { 288 INEXACT_MASK = 0x600, 289 UNDERFLOW_MASK = 0x010, 290 OVERFLOW_MASK = 0x008, 291 DIVBYZERO_MASK = 0x020, 292 INVALID_MASK = 0xF80 293 } 294 } else { // SPARC FSR is a 32bit register 295 //(64 bits for Sparc 7 & 8, but high 32 bits are uninteresting). 296 enum : int { 297 INEXACT_MASK = 0x020, 298 UNDERFLOW_MASK = 0x080, 299 OVERFLOW_MASK = 0x100, 300 DIVBYZERO_MASK = 0x040, 301 INVALID_MASK = 0x200 302 } 303 } 304 private: 305 @property static IeeeFlags getIeeeFlags() 306 { 307 // This is a highly time-critical operation, and 308 // should really be an intrinsic. 309 version(D_InlineAsm_X86) { 310 version(DMDWindows) { 311 // In this case, we 312 // take advantage of the fact that for DMD-Windows 313 // a struct containing only a int is returned in EAX. 314 asm { 315 fstsw AX; 316 // NOTE: If compiler supports SSE2, need to OR the result with 317 // the SSE2 status register. 318 // Clear all irrelevant bits 319 and EAX, 0x03D; 320 } 321 } 322 else { 323 IeeeFlags tmp1; 324 asm { 325 fstsw AX; 326 // NOTE: If compiler supports SSE2, need to OR the result with 327 // the SSE2 status register. 328 // Clear all irrelevant bits 329 and EAX, 0x03D; 330 mov tmp1, EAX; 331 } 332 return tmp1; 333 } 334 } else version (PPC) { 335 assert(0, "Not yet supported"); 336 } else { 337 /* SPARC: 338 int retval; 339 asm { st %fsr, retval; } 340 return retval; 341 */ 342 assert(0, "Not yet supported"); 343 } 344 } 345 @property static void resetIeeeFlags() 346 { 347 version(D_InlineAsm_X86) { 348 asm { 349 fnclex; 350 } 351 } else { 352 /* SPARC: 353 int tmpval; 354 asm { st %fsr, tmpval; } 355 tmpval &=0xFFFF_FC00; 356 asm { ld tmpval, %fsr; } 357 */ 358 assert(0, "Not yet supported"); 359 } 360 } 361 public: 362 /// The result cannot be represented exactly, so rounding occured. 363 /// (example: x = sin(0.1); } 364 @property int inexact() { return m_flags & INEXACT_MASK; } 365 /// A zero was generated by underflow (example: x = real.min_normal*real.epsilon/2;) 366 @property int underflow() { return m_flags & UNDERFLOW_MASK; } 367 /// An infinity was generated by overflow (example: x = real.max*2;) 368 @property int overflow() { return m_flags & OVERFLOW_MASK; } 369 /// An infinity was generated by division by zero (example: x = 3/0.0; ) 370 @property int divByZero() { return m_flags & DIVBYZERO_MASK; } 371 /// A machine NaN was generated. (example: x = real.infinity * 0.0; ) 372 @property int invalid() { return m_flags & INVALID_MASK; } 373 } 374 375 /// Return a snapshot of the current state of the floating-point status flags. 376 @property IeeeFlags ieeeFlags() { return IeeeFlags.getIeeeFlags(); } 377 378 /// Set all of the floating-point status flags to false. 379 void resetIeeeFlags() { IeeeFlags.resetIeeeFlags(); } 380 381 /** IEEE rounding modes. 382 * The default mode is ROUNDTONEAREST. 383 */ 384 enum RoundingMode : short { 385 ROUNDTONEAREST = 0x0000, 386 ROUNDDOWN = 0x0400, 387 ROUNDUP = 0x0800, 388 ROUNDTOZERO = 0x0C00 389 }; 390 391 /** Change the rounding mode used for all floating-point operations. 392 * 393 * Returns the old rounding mode. 394 * 395 * When changing the rounding mode, it is almost always necessary to restore it 396 * at the end of the function. Typical usage: 397 --- 398 auto oldrounding = setIeeeRounding(RoundingMode.ROUNDDOWN); 399 scope (exit) setIeeeRounding(oldrounding); 400 --- 401 */ 402 RoundingMode setIeeeRounding(RoundingMode roundingmode) { 403 version(D_InlineAsm_X86) { 404 // TODO: For SSE/SSE2, do we also need to set the SSE rounding mode? 405 short cont; 406 asm { 407 fstcw cont; 408 mov CX, cont; 409 mov AX, cont; 410 and EAX, 0x0C00; // Form the return value 411 and CX, 0xF3FF; 412 or CX, roundingmode; 413 mov cont, CX; 414 fldcw cont; 415 } 416 } else { 417 assert(0, "Not yet supported"); 418 } 419 } 420 421 /** Get the IEEE rounding mode which is in use. 422 * 423 */ 424 RoundingMode getIeeeRounding() { 425 version(D_InlineAsm_X86) { 426 // TODO: For SSE/SSE2, do we also need to check the SSE rounding mode? 427 short cont; 428 asm { 429 mov EAX, 0x0C00; 430 fstcw cont; 431 and AX, cont; 432 } 433 } else { 434 assert(0, "Not yet supported"); 435 } 436 } 437 438 debug(UnitTest) { 439 version(D_InlineAsm_X86) { // Won't work for anything else yet 440 unittest { 441 real a = 3.5; 442 resetIeeeFlags(); 443 assert(!ieeeFlags.divByZero); 444 a /= 0.0L; 445 assert(ieeeFlags.divByZero); 446 assert(a == real.infinity); 447 a *= 0.0L; 448 assert(ieeeFlags.invalid); 449 assert(isNaN(a)); 450 a = real.max; 451 a *= 2; 452 assert(ieeeFlags.overflow); 453 a = real.min_normal * real.epsilon; 454 a /= 99; 455 assert(ieeeFlags.underflow); 456 assert(ieeeFlags.inexact); 457 458 int r = getIeeeRounding(); 459 assert(r == RoundingMode.ROUNDTONEAREST); 460 } 461 } 462 } 463 464 // Note: Itanium supports more precision options than this. SSE/SSE2 does not support any. 465 enum PrecisionControl : short { 466 PRECISION80 = 0x300, 467 PRECISION64 = 0x200, 468 PRECISION32 = 0x000 469 }; 470 471 /** Set the number of bits of precision used by 'real'. 472 * 473 * Returns: the old precision. 474 * This is not supported on all platforms. 475 */ 476 PrecisionControl reduceRealPrecision(PrecisionControl prec) { 477 version(D_InlineAsm_X86) { 478 short cont; 479 asm { 480 fstcw cont; 481 mov CX, cont; 482 mov AX, cont; 483 and EAX, 0x0300; // Form the return value 484 and CX, 0xFCFF; 485 or CX, prec; 486 mov cont, CX; 487 fldcw cont; 488 } 489 } else { 490 assert(0, "Not yet supported"); 491 } 492 } 493 494 /********************************************************************* 495 * Separate floating point value into significand and exponent. 496 * 497 * Returns: 498 * Calculate and return $(I x) and $(I exp) such that 499 * value =$(I x)*2$(SUP exp) and 500 * .5 $(LT)= |$(I x)| $(LT) 1.0 501 * 502 * $(I x) has same sign as value. 503 * 504 * $(TABLE_SV 505 * $(TR $(TH value) $(TH returns) $(TH exp)) 506 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0)) 507 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max)) 508 * $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min)) 509 * $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min)) 510 * ) 511 */ 512 real frexp(real value, out int exp) 513 { 514 ushort* vu = cast(ushort*)&value; 515 long* vl = cast(long*)&value; 516 uint ex; 517 alias floatTraits!(real) F; 518 519 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 520 static if (real.mant_dig == 64) { // real80 521 if (ex) { // If exponent is non-zero 522 if (ex == F.EXPMASK) { // infinity or NaN 523 if (*vl & 0x7FFF_FFFF_FFFF_FFFF) { // NaN 524 *vl |= 0xC000_0000_0000_0000; // convert $(NAN)S to $(NAN)Q 525 exp = int.min; 526 } else if (vu[F.EXPPOS_SHORT] & 0x8000) { // negative infinity 527 exp = int.min; 528 } else { // positive infinity 529 exp = int.max; 530 } 531 } else { 532 exp = ex - F.EXPBIAS; 533 vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE); 534 } 535 } else if (!*vl) { 536 // value is +-0.0 537 exp = 0; 538 } else { 539 // denormal 540 value *= F.RECIP_EPSILON; 541 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 542 exp = ex - F.EXPBIAS - 63; 543 vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE); 544 } 545 return value; 546 } else static if (real.mant_dig == 113) { // quadruple 547 if (ex) { // If exponent is non-zero 548 if (ex == F.EXPMASK) { // infinity or NaN 549 if (vl[MANTISSA_LSB] |( vl[MANTISSA_MSB]&0x0000_FFFF_FFFF_FFFF)) { // NaN 550 vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000; // convert $(NAN)S to $(NAN)Q 551 exp = int.min; 552 } else if (vu[F.EXPPOS_SHORT] & 0x8000) { // negative infinity 553 exp = int.min; 554 } else { // positive infinity 555 exp = int.max; 556 } 557 } else { 558 exp = ex - F.EXPBIAS; 559 vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE); 560 } 561 } else if ((vl[MANTISSA_LSB] |(vl[MANTISSA_MSB]&0x0000_FFFF_FFFF_FFFF))==0) { 562 // value is +-0.0 563 exp = 0; 564 } else { 565 // denormal 566 value *= F.RECIP_EPSILON; 567 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; 568 exp = ex - F.EXPBIAS - 113; 569 vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE); 570 } 571 return value; 572 } else static if (real.mant_dig==53) { // real is double 573 if (ex) { // If exponent is non-zero 574 if (ex == F.EXPMASK) { // infinity or NaN 575 if (*vl==0x7FF0_0000_0000_0000) { // positive infinity 576 exp = int.max; 577 } else if (*vl==0xFFF0_0000_0000_0000) { // negative infinity 578 exp = int.min; 579 } else { // NaN 580 *vl |= 0x0008_0000_0000_0000; // convert $(NAN)S to $(NAN)Q 581 exp = int.min; 582 } 583 } else { 584 exp = (ex - F.EXPBIAS) >>> 4; 585 vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FE0; 586 } 587 } else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF)) { 588 // value is +-0.0 589 exp = 0; 590 } else { 591 // denormal 592 ushort sgn; 593 sgn = (0x8000 & vu[F.EXPPOS_SHORT])| 0x3FE0; 594 *vl &= 0x7FFF_FFFF_FFFF_FFFF; 595 596 int i = -0x3FD+11; 597 do { 598 i--; 599 *vl <<= 1; 600 } while (*vl > 0); 601 exp = i; 602 vu[F.EXPPOS_SHORT] = sgn; 603 } 604 return value; 605 }else { //static if(real.mant_dig==106) // doubledouble 606 assert(0, "Unsupported"); 607 } 608 } 609 610 debug(UnitTest) { 611 612 unittest 613 { 614 static real[3][] vals = // x,frexp,exp 615 [ 616 [0.0, 0.0, 0], 617 [-0.0, -0.0, 0], 618 [1.0, .5, 1], 619 [-1.0, -.5, 1], 620 [2.0, .5, 2], 621 [double.min_normal/2.0, .5, -1022], 622 [real.infinity,real.infinity,int.max], 623 [-real.infinity,-real.infinity,int.min], 624 ]; 625 626 int i; 627 int eptr; 628 real v = frexp(NaN(0xABC), eptr); 629 assert(isIdentical(NaN(0xABC), v)); 630 assert(eptr ==int.min); 631 v = frexp(-NaN(0xABC), eptr); 632 assert(isIdentical(-NaN(0xABC), v)); 633 assert(eptr ==int.min); 634 635 for (i = 0; i < vals.length; i++) { 636 real x = vals[i][0]; 637 real e = vals[i][1]; 638 int exp = cast(int)vals[i][2]; 639 v = frexp(x, eptr); 640 // printf("frexp(%La) = %La, should be %La, eptr = %d, should be %d\n", x, v, e, eptr, exp); 641 assert(isIdentical(e, v)); 642 assert(exp == eptr); 643 644 } 645 static if (real.mant_dig == 64) { 646 static real[3][] extendedvals = [ // x,frexp,exp 647 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal 648 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063], 649 [real.min_normal, .5, -16381], 650 [real.min_normal/2.0L, .5, -16382] // denormal 651 ]; 652 653 for (i = 0; i < extendedvals.length; i++) { 654 real x = extendedvals[i][0]; 655 real e = extendedvals[i][1]; 656 int exp = cast(int)extendedvals[i][2]; 657 v = frexp(x, eptr); 658 assert(isIdentical(e, v)); 659 assert(exp == eptr); 660 661 } 662 } 663 } 664 } 665 666 /** 667 * Compute n * 2$(SUP exp) 668 * References: frexp 669 */ 670 real ldexp(real n, int exp) /* intrinsic */ 671 { 672 version(Naked_D_InlineAsm_X86) 673 { 674 asm { 675 fild exp; 676 fld n; 677 fscale; 678 fstp ST(1);//, ST(0); 679 } 680 } 681 else 682 { 683 return tango.stdc.math.ldexpl(n, exp); 684 } 685 } 686 687 /****************************************** 688 * Extracts the exponent of x as a signed integral value. 689 * 690 * If x is not a special value, the result is the same as 691 * $(D cast(int)logb(x)). 692 * 693 * Remarks: This function is consistent with IEEE754R, but it 694 * differs from the C function of the same name 695 * in the return value of infinity. (in C, ilogb(real.infinity)== int.max). 696 * Note that the special return values may all be equal. 697 * 698 * $(TABLE_SV 699 * $(TR $(TH x) $(TH ilogb(x)) $(TH Invalid?)) 700 * $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes)) 701 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD FP_ILOGBINFINITY) $(TD yes)) 702 * $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD yes)) 703 * ) 704 */ 705 int ilogb(real x) 706 { 707 version(Naked_D_InlineAsm_X86) 708 { 709 int y; 710 asm { 711 fld x; 712 fxtract; 713 fstp ST(0); // drop significand 714 fistp y; // and return the exponent 715 } 716 return y; 717 } else static if (real.mant_dig==64) { // 80-bit reals 718 alias floatTraits!(real) F; 719 short e = cast(short)((cast(short *)&x)[F.EXPPOS_SHORT] & F.EXPMASK); 720 if (e == F.EXPMASK) { 721 // BUG: should also set the invalid exception 722 ulong s = *cast(ulong *)&x; 723 if (s == 0x8000_0000_0000_0000) { 724 return FP_ILOGBINFINITY; 725 } 726 else return FP_ILOGBNAN; 727 } 728 if (e==0) { 729 ulong s = *cast(ulong *)&x; 730 if (s == 0x0000_0000_0000_0000) { 731 // BUG: should also set the invalid exception 732 return FP_ILOGB0; 733 } 734 // Denormals 735 x *= F.RECIP_EPSILON; 736 short f = (cast(short *)&x)[F.EXPPOS_SHORT]; 737 return -0x3FFF - (63-f); 738 } 739 return e - 0x3FFF; 740 } else { 741 return tango.stdc.math.ilogbl(x); 742 } 743 } 744 745 version (X86) 746 { 747 enum int FP_ILOGB0 = -int.max-1; 748 enum int FP_ILOGBNAN = -int.max-1; 749 enum int FP_ILOGBINFINITY = -int.max-1; 750 } else { 751 alias tango.stdc.math.FP_ILOGB0 FP_ILOGB0; 752 alias tango.stdc.math.FP_ILOGBNAN FP_ILOGBNAN; 753 enum int FP_ILOGBINFINITY = int.max; 754 } 755 756 debug(UnitTest) { 757 unittest { 758 assert(ilogb(1.0) == 0); 759 assert(ilogb(65536) == 16); 760 assert(ilogb(-65536) == 16); 761 assert(ilogb(1.0 / 65536) == -16); 762 assert(ilogb(real.nan) == FP_ILOGBNAN); 763 assert(ilogb(0.0) == FP_ILOGB0); 764 assert(ilogb(-0.0) == FP_ILOGB0); 765 // denormal 766 assert(ilogb(0.125 * real.min_normal) == real.min_exp - 4); 767 assert(ilogb(real.infinity) == FP_ILOGBINFINITY); 768 } 769 } 770 771 /***************************************** 772 * Extracts the exponent of x as a signed integral value. 773 * 774 * If x is subnormal, it is treated as if it were normalized. 775 * For a positive, finite x: 776 * 777 * 1 $(LT)= $(I x) * FLT_RADIX$(SUP -logb(x)) $(LT) FLT_RADIX 778 * 779 * $(TABLE_SV 780 * $(TR $(TH x) $(TH logb(x)) $(TH divide by 0?) ) 781 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no)) 782 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) ) 783 * ) 784 */ 785 real logb(real x) 786 { 787 version(Naked_D_InlineAsm_X86) 788 { 789 asm { 790 fld x; 791 fxtract; 792 fstp ST(0);//, ST; // drop significand 793 } 794 } else { 795 return tango.stdc.math.logbl(x); 796 } 797 } 798 799 debug(UnitTest) { 800 unittest { 801 assert(logb(real.infinity)== real.infinity); 802 assert(isIdentical(logb(NaN(0xFCD)), NaN(0xFCD))); 803 assert(logb(1.0)== 0.0); 804 assert(logb(-65536) == 16); 805 assert(logb(0.0)== -real.infinity); 806 assert(ilogb(0.125*real.min_normal) == real.min_exp-4); 807 } 808 } 809 810 /************************************* 811 * Efficiently calculates x * 2$(SUP n). 812 * 813 * scalbn handles underflow and overflow in 814 * the same fashion as the basic arithmetic operators. 815 * 816 * $(TABLE_SV 817 * $(TR $(TH x) $(TH scalb(x))) 818 * $(TR $(TD $(PLUSMNINF)) $(TD $(PLUSMNINF)) ) 819 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) ) 820 * ) 821 */ 822 real scalbn(real x, int n) 823 { 824 version(Naked_D_InlineAsm_X86) 825 { 826 asm { 827 fild n; 828 fld x; 829 fscale; 830 fstp ST(1);//, ST; 831 } 832 } else { 833 // NOTE: Not implemented in DMD 834 return tango.stdc.math.scalbnl(x, n); 835 } 836 } 837 838 debug(UnitTest) { 839 unittest { 840 assert(scalbn(-real.infinity, 5) == -real.infinity); 841 assert(isIdentical(scalbn(NaN(0xABC),7), NaN(0xABC))); 842 } 843 } 844 845 /** 846 * Returns the positive difference between x and y. 847 * 848 * If either of x or y is $(NAN), it will be returned. 849 * Returns: 850 * $(TABLE_SV 851 * $(SVH Arguments, fdim(x, y)) 852 * $(SV x $(GT) y, x - y) 853 * $(SV x $(LT)= y, +0.0) 854 * ) 855 */ 856 real fdim(real x, real y) 857 { 858 return (x.isNaN || y.isNaN || x > y) ? x - y : +0.0; 859 } 860 861 debug(UnitTest) { 862 unittest { 863 assert(isIdentical(fdim(NaN(0xABC), 58.2), NaN(0xABC))); 864 } 865 } 866 867 /******************************* 868 * Returns |x| 869 * 870 * $(TABLE_SV 871 * $(TR $(TH x) $(TH fabs(x))) 872 * $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) ) 873 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) ) 874 * ) 875 */ 876 real fabs(real x) /* intrinsic */ 877 { 878 version(D_InlineAsm_X86) 879 { 880 asm { 881 fld x; 882 fabs; 883 } 884 } 885 else 886 { 887 return tango.stdc.math.fabsl(x); 888 } 889 } 890 891 unittest { 892 assert(isIdentical(fabs(NaN(0xABC)), NaN(0xABC))); 893 } 894 895 /** 896 * Returns (x * y) + z, rounding only once according to the 897 * current rounding mode. 898 * 899 * BUGS: Not currently implemented - rounds twice. 900 */ 901 real fma(float x, float y, float z) 902 { 903 return (x * y) + z; 904 } 905 906 /** 907 * Calculate cos(y) + i sin(y). 908 * 909 * On x86 CPUs, this is a very efficient operation; 910 * almost twice as fast as calculating sin(y) and cos(y) 911 * seperately, and is the preferred method when both are required. 912 */ 913 creal expi(real y) 914 { 915 version(Naked_D_InlineAsm_X86) 916 { 917 asm { 918 fld y; 919 fsincos; 920 fxch ST(1), ST(0); 921 } 922 } 923 else 924 { 925 return tango.stdc.math.cosl(y) + tango.stdc.math.sinl(y)*1i; 926 } 927 } 928 929 debug(UnitTest) { 930 unittest 931 { 932 assert(expi(1.3e5L) == tango.stdc.math.cosl(1.3e5L) + tango.stdc.math.sinl(1.3e5L) * 1i); 933 assert(expi(0.0L) == 1L + 0.0Li); 934 } 935 } 936 937 /********************************* 938 * Returns !=0 if e is a NaN. 939 */ 940 941 int isNaN(real x) 942 { 943 alias floatTraits!(real) F; 944 static if (real.mant_dig==53) { // double 945 ulong* p = cast(ulong *)&x; 946 return ((*p & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000) && *p & 0x000F_FFFF_FFFF_FFFF; 947 } else static if (real.mant_dig==64) { // real80 948 ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]; 949 ulong* ps = cast(ulong *)&x; 950 return e == F.EXPMASK && 951 *ps & 0x7FFF_FFFF_FFFF_FFFF; // not infinity 952 } else static if (real.mant_dig==113) { // quadruple 953 ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]; 954 ulong* ps = cast(ulong *)&x; 955 return e == F.EXPMASK && 956 (ps[MANTISSA_LSB] | (ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF))!=0; 957 } else { 958 return x!=x; 959 } 960 } 961 962 963 debug(UnitTest) { 964 unittest 965 { 966 assert(isNaN(float.nan)); 967 assert(isNaN(-double.nan)); 968 assert(isNaN(real.nan)); 969 970 assert(!isNaN(53.6)); 971 assert(!isNaN(float.infinity)); 972 } 973 } 974 975 /** 976 * Returns !=0 if x is normalized. 977 * 978 * (Need one for each format because subnormal 979 * floats might be converted to normal reals) 980 */ 981 int isNormal(X)(X x) 982 { 983 alias floatTraits!(X) F; 984 985 static if(real.mant_dig==106) { // doubledouble 986 // doubledouble is normal if the least significant part is normal. 987 return isNormal((cast(double*)&x)[MANTISSA_LSB]); 988 } else { 989 ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]; 990 return (e != F.EXPMASK && e!=0); 991 } 992 } 993 994 debug(UnitTest) { 995 unittest 996 { 997 float f = 3; 998 double d = 500; 999 real e = 10e+48; 1000 1001 assert(isNormal(f)); 1002 assert(isNormal(d)); 1003 assert(isNormal(e)); 1004 f=d=e=0; 1005 assert(!isNormal(f)); 1006 assert(!isNormal(d)); 1007 assert(!isNormal(e)); 1008 assert(!isNormal(real.infinity)); 1009 assert(isNormal(-real.max)); 1010 assert(!isNormal(real.min_normal/4)); 1011 1012 } 1013 } 1014 1015 /********************************* 1016 * Is the binary representation of x identical to y? 1017 * 1018 * Same as ==, except that positive and negative zero are not identical, 1019 * and two $(NAN)s are identical if they have the same 'payload'. 1020 */ 1021 1022 bool isIdentical(real x, real y) 1023 { 1024 // We're doing a bitwise comparison so the endianness is irrelevant. 1025 long* pxs = cast(long *)&x; 1026 long* pys = cast(long *)&y; 1027 static if (real.mant_dig == 53){ //double 1028 return pxs[0] == pys[0]; 1029 } else static if (real.mant_dig == 113 || real.mant_dig==106) { 1030 // quadruple or doubledouble 1031 return pxs[0] == pys[0] && pxs[1] == pys[1]; 1032 } else { // real80 1033 ushort* pxe = cast(ushort *)&x; 1034 ushort* pye = cast(ushort *)&y; 1035 return pxe[4] == pye[4] && pxs[0] == pys[0]; 1036 } 1037 } 1038 1039 /** ditto */ 1040 bool isIdentical(ireal x, ireal y) { 1041 return isIdentical(x.im, y.im); 1042 } 1043 1044 /** ditto */ 1045 bool isIdentical(creal x, creal y) { 1046 return isIdentical(x.re, y.re) && isIdentical(x.im, y.im); 1047 } 1048 1049 debug(UnitTest) { 1050 unittest { 1051 assert(isIdentical(0.0, 0.0)); 1052 assert(!isIdentical(0.0, -0.0)); 1053 assert(isIdentical(NaN(0xABC), NaN(0xABC))); 1054 assert(!isIdentical(NaN(0xABC), NaN(218))); 1055 assert(isIdentical(1.234e56, 1.234e56)); 1056 assert(isNaN(NaN(0x12345))); 1057 assert(isIdentical(3.1 + NaN(0xDEF) * 1i, 3.1 + NaN(0xDEF)*1i)); 1058 assert(!isIdentical(3.1+0.0i, 3.1-0i)); 1059 assert(!isIdentical(0.0i, 2.5e58i)); 1060 } 1061 } 1062 1063 /********************************* 1064 * Is number subnormal? (Also called "denormal".) 1065 * Subnormals have a 0 exponent and a 0 most significant significand bit, 1066 * but are non-zero. 1067 */ 1068 1069 /* Need one for each format because subnormal floats might 1070 * be converted to normal reals. 1071 */ 1072 1073 int isSubnormal(float f) 1074 { 1075 uint *p = cast(uint *)&f; 1076 return (*p & 0x7F80_0000) == 0 && *p & 0x007F_FFFF; 1077 } 1078 1079 debug(UnitTest) { 1080 unittest 1081 { 1082 float f = -float.min_normal; 1083 assert(!isSubnormal(f)); 1084 f/=4; 1085 assert(isSubnormal(f)); 1086 } 1087 } 1088 1089 /// ditto 1090 1091 int isSubnormal(double d) 1092 { 1093 uint *p = cast(uint *)&d; 1094 return (p[MANTISSA_MSB] & 0x7FF0_0000) == 0 && (p[MANTISSA_LSB] || p[MANTISSA_MSB] & 0x000F_FFFF); 1095 } 1096 1097 debug(UnitTest) { 1098 unittest 1099 { 1100 double f; 1101 1102 for (f = 1; !isSubnormal(f); f /= 2) 1103 assert(f != 0); 1104 } 1105 } 1106 1107 /// ditto 1108 1109 int isSubnormal(real x) 1110 { 1111 alias floatTraits!(real) F; 1112 static if (real.mant_dig == 53) { // double 1113 return isSubnormal(cast(double)x); 1114 } else static if (real.mant_dig == 113) { // quadruple 1115 ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]; 1116 long* ps = cast(long *)&x; 1117 return (e == 0 && (((ps[MANTISSA_LSB]|(ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF))) !=0)); 1118 } else static if (real.mant_dig==64) { // real80 1119 ushort* pe = cast(ushort *)&x; 1120 long* ps = cast(long *)&x; 1121 1122 return (pe[F.EXPPOS_SHORT] & F.EXPMASK) == 0 && *ps > 0; 1123 } else { // double double 1124 return isSubnormal((cast(double*)&x)[MANTISSA_MSB]); 1125 } 1126 } 1127 1128 debug(UnitTest) { 1129 unittest 1130 { 1131 real f; 1132 1133 for (f = 1; !isSubnormal(f); f /= 2) 1134 assert(f != 0); 1135 } 1136 } 1137 1138 /********************************* 1139 * Return !=0 if x is $(PLUSMN)0. 1140 * 1141 * Does not affect any floating-point flags 1142 */ 1143 int isZero(real x) 1144 { 1145 alias floatTraits!(real) F; 1146 static if (real.mant_dig == 53) { // double 1147 return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF) == 0; 1148 } else static if (real.mant_dig == 113) { // quadruple 1149 long* ps = cast(long *)&x; 1150 return (ps[MANTISSA_LSB] | (ps[MANTISSA_MSB]& 0x7FFF_FFFF_FFFF_FFFF)) == 0; 1151 } else { // real80 1152 ushort* pe = cast(ushort *)&x; 1153 ulong* ps = cast(ulong *)&x; 1154 return (pe[F.EXPPOS_SHORT] & F.EXPMASK) == 0 && *ps == 0; 1155 } 1156 } 1157 1158 debug(UnitTest) { 1159 unittest 1160 { 1161 assert(isZero(0.0)); 1162 assert(isZero(-0.0)); 1163 assert(!isZero(2.5)); 1164 assert(!isZero(real.min_normal / 1000)); 1165 } 1166 } 1167 1168 /********************************* 1169 * Return !=0 if e is $(PLUSMNINF);. 1170 */ 1171 1172 int isInfinity(real x) 1173 { 1174 alias floatTraits!(real) F; 1175 static if (real.mant_dig == 53) { // double 1176 return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF8_0000_0000_0000; 1177 } else static if(real.mant_dig == 106) { //doubledouble 1178 return (((cast(ulong *)&x)[MANTISSA_MSB]) & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF8_0000_0000_0000; 1179 } else static if (real.mant_dig == 113) { // quadruple 1180 long* ps = cast(long *)&x; 1181 return (ps[MANTISSA_LSB] == 0) 1182 && (ps[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_0000_0000_0000; 1183 } else { // real80 1184 ushort e = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]); 1185 ulong* ps = cast(ulong *)&x; 1186 1187 return e == F.EXPMASK && *ps == 0x8000_0000_0000_0000; 1188 } 1189 } 1190 1191 debug(UnitTest) { 1192 unittest 1193 { 1194 assert(isInfinity(float.infinity)); 1195 assert(!isInfinity(float.nan)); 1196 assert(isInfinity(double.infinity)); 1197 assert(isInfinity(-real.infinity)); 1198 1199 assert(isInfinity(-1.0 / 0.0)); 1200 } 1201 } 1202 1203 /** 1204 * Calculate the next largest floating point value after x. 1205 * 1206 * Return the least number greater than x that is representable as a real; 1207 * thus, it gives the next point on the IEEE number line. 1208 * 1209 * $(TABLE_SV 1210 * $(SVH x, nextUp(x) ) 1211 * $(SV -$(INFIN), -real.max ) 1212 * $(SV $(PLUSMN)0.0, real.min_normal*real.epsilon ) 1213 * $(SV real.max, $(INFIN) ) 1214 * $(SV $(INFIN), $(INFIN) ) 1215 * $(SV $(NAN), $(NAN) ) 1216 * ) 1217 * 1218 * Remarks: 1219 * This function is included in the IEEE 754-2008 standard. 1220 * 1221 * nextDoubleUp and nextFloatUp are the corresponding functions for 1222 * the IEEE double and IEEE float number lines. 1223 */ 1224 real nextUp(real x) 1225 { 1226 alias floatTraits!(real) F; 1227 static if (real.mant_dig == 53) { // double 1228 return nextDoubleUp(x); 1229 } else static if(real.mant_dig==113) { // quadruple 1230 ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]; 1231 if (e == F.EXPMASK) { // NaN or Infinity 1232 if (x == -real.infinity) return -real.max; 1233 return x; // +Inf and NaN are unchanged. 1234 } 1235 ulong* ps = cast(ulong *)&e; 1236 if (ps[MANTISSA_LSB] & 0x8000_0000_0000_0000) { // Negative number 1237 if (ps[MANTISSA_LSB]==0 && ps[MANTISSA_MSB] == 0x8000_0000_0000_0000) { // it was negative zero 1238 ps[MANTISSA_LSB] = 0x0000_0000_0000_0001; // change to smallest subnormal 1239 ps[MANTISSA_MSB] = 0; 1240 return x; 1241 } 1242 --*ps; 1243 if (ps[MANTISSA_LSB]==0) --ps[MANTISSA_MSB]; 1244 } else { // Positive number 1245 ++ps[MANTISSA_LSB]; 1246 if (ps[MANTISSA_LSB]==0) ++ps[MANTISSA_MSB]; 1247 } 1248 return x; 1249 1250 } else static if(real.mant_dig==64){ // real80 1251 // For 80-bit reals, the "implied bit" is a nuisance... 1252 ushort *pe = cast(ushort *)&x; 1253 ulong *ps = cast(ulong *)&x; 1254 1255 if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK) { 1256 // First, deal with NANs and infinity 1257 if (x == -real.infinity) return -real.max; 1258 return x; // +Inf and NaN are unchanged. 1259 } 1260 if (pe[F.EXPPOS_SHORT] & 0x8000) { // Negative number -- need to decrease the significand 1261 --*ps; 1262 // Need to mask with 0x7FFF... so subnormals are treated correctly. 1263 if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF) { 1264 if (pe[F.EXPPOS_SHORT] == 0x8000) { // it was negative zero 1265 *ps = 1; 1266 pe[F.EXPPOS_SHORT] = 0; // smallest subnormal. 1267 return x; 1268 } 1269 --pe[F.EXPPOS_SHORT]; 1270 if (pe[F.EXPPOS_SHORT] == 0x8000) { 1271 return x; // it's become a subnormal, implied bit stays low. 1272 } 1273 *ps = 0xFFFF_FFFF_FFFF_FFFF; // set the implied bit 1274 return x; 1275 } 1276 return x; 1277 } else { 1278 // Positive number -- need to increase the significand. 1279 // Works automatically for positive zero. 1280 ++*ps; 1281 if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0) { 1282 // change in exponent 1283 ++pe[F.EXPPOS_SHORT]; 1284 *ps = 0x8000_0000_0000_0000; // set the high bit 1285 } 1286 } 1287 return x; 1288 } else { // doubledouble 1289 assert(0, "Not implemented"); 1290 } 1291 } 1292 1293 /** ditto */ 1294 double nextDoubleUp(double x) 1295 { 1296 ulong *ps = cast(ulong *)&x; 1297 1298 if ((*ps & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000) { 1299 // First, deal with NANs and infinity 1300 if (x == -x.infinity) return -x.max; 1301 return x; // +INF and NAN are unchanged. 1302 } 1303 if (*ps & 0x8000_0000_0000_0000) { // Negative number 1304 if (*ps == 0x8000_0000_0000_0000) { // it was negative zero 1305 *ps = 0x0000_0000_0000_0001; // change to smallest subnormal 1306 return x; 1307 } 1308 --*ps; 1309 } else { // Positive number 1310 ++*ps; 1311 } 1312 return x; 1313 } 1314 1315 /** ditto */ 1316 float nextFloatUp(float x) 1317 { 1318 uint *ps = cast(uint *)&x; 1319 1320 if ((*ps & 0x7F80_0000) == 0x7F80_0000) { 1321 // First, deal with NANs and infinity 1322 if (x == -x.infinity) return -x.max; 1323 return x; // +INF and NAN are unchanged. 1324 } 1325 if (*ps & 0x8000_0000) { // Negative number 1326 if (*ps == 0x8000_0000) { // it was negative zero 1327 *ps = 0x0000_0001; // change to smallest subnormal 1328 return x; 1329 } 1330 --*ps; 1331 } else { // Positive number 1332 ++*ps; 1333 } 1334 return x; 1335 } 1336 1337 debug(UnitTest) { 1338 unittest { 1339 static if (real.mant_dig == 64) { 1340 1341 // Tests for 80-bit reals 1342 1343 assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC))); 1344 // negative numbers 1345 assert( nextUp(-real.infinity) == -real.max ); 1346 assert( nextUp(-1-real.epsilon) == -1.0 ); 1347 assert( nextUp(-2) == -2.0 + real.epsilon); 1348 // denormals and zero 1349 assert( nextUp(-real.min_normal) == -real.min_normal*(1-real.epsilon) ); 1350 assert( nextUp(-real.min_normal*(1-real.epsilon) == -real.min_normal*(1-2*real.epsilon)) ); 1351 assert( isIdentical(-0.0L, nextUp(-real.min_normal*real.epsilon)) ); 1352 assert( nextUp(-0.0) == real.min_normal*real.epsilon ); 1353 assert( nextUp(0.0) == real.min_normal*real.epsilon ); 1354 assert( nextUp(real.min_normal*(1-real.epsilon)) == real.min_normal ); 1355 assert( nextUp(real.min_normal) == real.min_normal*(1+real.epsilon) ); 1356 // positive numbers 1357 assert( nextUp(1) == 1.0 + real.epsilon ); 1358 assert( nextUp(2.0-real.epsilon) == 2.0 ); 1359 assert( nextUp(real.max) == real.infinity ); 1360 assert( nextUp(real.infinity)==real.infinity ); 1361 } 1362 1363 assert(isIdentical(nextDoubleUp(NaN(0xABC)), NaN(0xABC))); 1364 // negative numbers 1365 assert( nextDoubleUp(-double.infinity) == -double.max ); 1366 assert( nextDoubleUp(-1-double.epsilon) == -1.0 ); 1367 assert( nextDoubleUp(-2) == -2.0 + double.epsilon); 1368 // denormals and zero 1369 1370 assert( nextDoubleUp(-double.min_normal) == -double.min_normal*(1-double.epsilon) ); 1371 assert( nextDoubleUp(-double.min_normal*(1-double.epsilon) == -double.min_normal*(1-2*double.epsilon)) ); 1372 assert( isIdentical(-0.0, nextDoubleUp(-double.min_normal*double.epsilon)) ); 1373 assert( nextDoubleUp(0.0) == double.min_normal*double.epsilon ); 1374 assert( nextDoubleUp(-0.0) == double.min_normal*double.epsilon ); 1375 assert( nextDoubleUp(double.min_normal*(1-double.epsilon)) == double.min_normal ); 1376 assert( nextDoubleUp(double.min_normal) == double.min_normal*(1+double.epsilon) ); 1377 // positive numbers 1378 assert( nextDoubleUp(1) == 1.0 + double.epsilon ); 1379 assert( nextDoubleUp(2.0-double.epsilon) == 2.0 ); 1380 assert( nextDoubleUp(double.max) == double.infinity ); 1381 1382 assert(isIdentical(nextFloatUp(NaN(0xABC)), NaN(0xABC))); 1383 assert( nextFloatUp(-float.min_normal) == -float.min_normal*(1-float.epsilon) ); 1384 assert( nextFloatUp(1.0) == 1.0+float.epsilon ); 1385 assert( nextFloatUp(-0.0) == float.min_normal*float.epsilon); 1386 assert( nextFloatUp(float.infinity)==float.infinity ); 1387 1388 assert(nextDown(1.0+real.epsilon)==1.0); 1389 assert(nextDoubleDown(1.0+double.epsilon)==1.0); 1390 assert(nextFloatDown(1.0+float.epsilon)==1.0); 1391 assert(nextafter(1.0+real.epsilon, -real.infinity)==1.0); 1392 } 1393 } 1394 1395 package { 1396 /** Reduces the magnitude of x, so the bits in the lower half of its significand 1397 * are all zero. Returns the amount which needs to be added to x to restore its 1398 * initial value; this amount will also have zeros in all bits in the lower half 1399 * of its significand. 1400 */ 1401 X splitSignificand(X)(ref X x) 1402 { 1403 if (isNaN(x) || isInfinity(x)) return 0; // don't change NaN or infinity 1404 X y = x; // copy the original value 1405 static if (X.mant_dig == float.mant_dig) { 1406 uint *ps = cast(uint *)&x; 1407 (*ps) &= 0xFFFF_FC00; 1408 } else static if (X.mant_dig == 53) { 1409 ulong *ps = cast(ulong *)&x; 1410 (*ps) &= 0xFFFF_FFFF_FC00_0000L; 1411 } else static if (X.mant_dig == 64){ // 80-bit real 1412 // An x87 real80 has 63 bits, because the 'implied' bit is stored explicitly. 1413 // This is annoying, because it means the significand cannot be 1414 // precisely halved. Instead, we split it into 31+32 bits. 1415 ulong *ps = cast(ulong *)&x; 1416 (*ps) &= 0xFFFF_FFFF_0000_0000L; 1417 } else static if (X.mant_dig==113) { // quadruple 1418 ulong *ps = cast(ulong *)&x; 1419 ps[MANTISSA_LSB] &= 0xFF00_0000_0000_0000L; 1420 } 1421 //else static assert(0, "Unsupported size"); 1422 1423 return y - x; 1424 } 1425 1426 unittest { 1427 double x = -0x1.234_567A_AAAA_AAp+250; 1428 double y = splitSignificand(x); 1429 assert(x == -0x1.234_5678p+250); 1430 assert(y == -0x0.000_000A_AAAA_A8p+248); 1431 assert(x + y == -0x1.234_567A_AAAA_AAp+250); 1432 } 1433 } 1434 1435 /** 1436 * Calculate the next smallest floating point value before x. 1437 * 1438 * Return the greatest number less than x that is representable as a real; 1439 * thus, it gives the previous point on the IEEE number line. 1440 * 1441 * $(TABLE_SV 1442 * $(SVH x, nextDown(x) ) 1443 * $(SV $(INFIN), real.max ) 1444 * $(SV $(PLUSMN)0.0, -real.min_normal*real.epsilon ) 1445 * $(SV -real.max, -$(INFIN) ) 1446 * $(SV -$(INFIN), -$(INFIN) ) 1447 * $(SV $(NAN), $(NAN) ) 1448 * ) 1449 * 1450 * Remarks: 1451 * This function is included in the IEEE 754-2008 standard. 1452 * 1453 * nextDoubleDown and nextFloatDown are the corresponding functions for 1454 * the IEEE double and IEEE float number lines. 1455 */ 1456 real nextDown(real x) 1457 { 1458 return -nextUp(-x); 1459 } 1460 1461 /** ditto */ 1462 double nextDoubleDown(double x) 1463 { 1464 return -nextDoubleUp(-x); 1465 } 1466 1467 /** ditto */ 1468 float nextFloatDown(float x) 1469 { 1470 return -nextFloatUp(-x); 1471 } 1472 1473 debug(UnitTest) { 1474 unittest { 1475 assert( nextDown(1.0 + real.epsilon) == 1.0); 1476 } 1477 } 1478 1479 /** 1480 * Calculates the next representable value after x in the direction of y. 1481 * 1482 * If y > x, the result will be the next largest floating-point value; 1483 * if y < x, the result will be the next smallest value. 1484 * If x == y, the result is y. 1485 * 1486 * Remarks: 1487 * This function is not generally very useful; it's almost always better to use 1488 * the faster functions nextUp() or nextDown() instead. 1489 * 1490 * IEEE 754 requirements not implemented: 1491 * The FE_INEXACT and FE_OVERFLOW exceptions will be raised if x is finite and 1492 * the function result is infinite. The FE_INEXACT and FE_UNDERFLOW 1493 * exceptions will be raised if the function value is subnormal, and x is 1494 * not equal to y. 1495 */ 1496 real nextafter(real x, real y) 1497 { 1498 if (x==y) return y; 1499 return (y>x) ? nextUp(x) : nextDown(x); 1500 } 1501 1502 /************************************** 1503 * To what precision is x equal to y? 1504 * 1505 * Returns: the number of significand bits which are equal in x and y. 1506 * eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision. 1507 * 1508 * $(TABLE_SV 1509 * $(SVH3 x, y, feqrel(x, y) ) 1510 * $(SV3 x, x, typeof(x).mant_dig ) 1511 * $(SV3 x, $(GT)= 2*x, 0 ) 1512 * $(SV3 x, $(LE)= x/2, 0 ) 1513 * $(SV3 $(NAN), any, 0 ) 1514 * $(SV3 any, $(NAN), 0 ) 1515 * ) 1516 * 1517 * Remarks: 1518 * This is a very fast operation, suitable for use in speed-critical code. 1519 */ 1520 int feqrel(X)(X x, X y) 1521 { 1522 /* Public Domain. Author: Don Clugston, 18 Aug 2005. 1523 */ 1524 static assert(is(X==real) || is(X==double) || is(X==float), "Only float, double, and real are supported by feqrel"); 1525 1526 static if (X.mant_dig == 106) { // doubledouble. 1527 int a = feqrel(cast(double*)(&x)[MANTISSA_MSB], cast(double*)(&y)[MANTISSA_MSB]); 1528 if (a != double.mant_dig) return a; 1529 return double.mant_dig + feqrel(cast(double*)(&x)[MANTISSA_LSB], cast(double*)(&y)[MANTISSA_LSB]); 1530 } else static if (X.mant_dig==64 || X.mant_dig==113 1531 || X.mant_dig==53 || X.mant_dig == 24) { 1532 if (x == y) return X.mant_dig; // ensure diff!=0, cope with INF. 1533 1534 X diff = fabs(x - y); 1535 1536 ushort *pa = cast(ushort *)(&x); 1537 ushort *pb = cast(ushort *)(&y); 1538 ushort *pd = cast(ushort *)(&diff); 1539 1540 alias floatTraits!(X) F; 1541 1542 // The difference in abs(exponent) between x or y and abs(x-y) 1543 // is equal to the number of significand bits of x which are 1544 // equal to y. If negative, x and y have different exponents. 1545 // If positive, x and y are equal to 'bitsdiff' bits. 1546 // AND with 0x7FFF to form the absolute value. 1547 // To avoid out-by-1 errors, we subtract 1 so it rounds down 1548 // if the exponents were different. This means 'bitsdiff' is 1549 // always 1 lower than we want, except that if bitsdiff==0, 1550 // they could have 0 or 1 bits in common. 1551 1552 static if (X.mant_dig==64 || X.mant_dig==113) { // real80 or quadruple 1553 int bitsdiff = ( ((pa[F.EXPPOS_SHORT] & F.EXPMASK) 1554 + (pb[F.EXPPOS_SHORT]& F.EXPMASK) 1555 - (0x8000-F.EXPMASK))>>1) 1556 - pd[F.EXPPOS_SHORT]; 1557 } else static if (X.mant_dig==53) { // double 1558 int bitsdiff = (( ((pa[F.EXPPOS_SHORT] & F.EXPMASK) 1559 + (pb[F.EXPPOS_SHORT] & F.EXPMASK) 1560 - (0x8000-F.EXPMASK))>>1) 1561 - (pd[F.EXPPOS_SHORT] & F.EXPMASK))>>4; 1562 } else static if (X.mant_dig == 24) { // float 1563 int bitsdiff = (( ((pa[F.EXPPOS_SHORT] & F.EXPMASK) 1564 + (pb[F.EXPPOS_SHORT] & F.EXPMASK) 1565 - (0x8000-F.EXPMASK))>>1) 1566 - (pd[F.EXPPOS_SHORT] & F.EXPMASK))>>7; 1567 } 1568 if (pd[F.EXPPOS_SHORT] == 0) 1569 { // Difference is denormal 1570 // For denormals, we need to add the number of zeros that 1571 // lie at the start of diff's significand. 1572 // We do this by multiplying by 2^real.mant_dig 1573 diff *= F.RECIP_EPSILON; 1574 return bitsdiff + X.mant_dig - pd[F.EXPPOS_SHORT]; 1575 } 1576 1577 if (bitsdiff > 0) 1578 return bitsdiff + 1; // add the 1 we subtracted before 1579 1580 // Avoid out-by-1 errors when factor is almost 2. 1581 static if (X.mant_dig==64 || X.mant_dig==113) { // real80 or quadruple 1582 return (bitsdiff == 0) ? (pa[F.EXPPOS_SHORT] == pb[F.EXPPOS_SHORT]) : 0; 1583 } else static if (X.mant_dig == 53 || X.mant_dig == 24) { // double or float 1584 return (bitsdiff == 0 && !((pa[F.EXPPOS_SHORT] ^ pb[F.EXPPOS_SHORT])& F.EXPMASK)) ? 1 : 0; 1585 } 1586 } else { 1587 assert(0, "Unsupported"); 1588 } 1589 } 1590 1591 debug(UnitTest) { 1592 unittest 1593 { 1594 // Exact equality 1595 assert(feqrel(real.max,real.max)==real.mant_dig); 1596 assert(feqrel(0.0L,0.0L)==real.mant_dig); 1597 assert(feqrel(7.1824L,7.1824L)==real.mant_dig); 1598 assert(feqrel(real.infinity,real.infinity)==real.mant_dig); 1599 1600 // a few bits away from exact equality 1601 real w=1; 1602 for (int i=1; i<real.mant_dig-1; ++i) { 1603 assert(feqrel(1+w*real.epsilon,1.0L)==real.mant_dig-i); 1604 assert(feqrel(1-w*real.epsilon,1.0L)==real.mant_dig-i); 1605 assert(feqrel(1.0L,1+(w-1)*real.epsilon)==real.mant_dig-i+1); 1606 w*=2; 1607 } 1608 assert(feqrel(1.5+real.epsilon,1.5L)==real.mant_dig-1); 1609 assert(feqrel(1.5-real.epsilon,1.5L)==real.mant_dig-1); 1610 assert(feqrel(1.5-real.epsilon,1.5+real.epsilon)==real.mant_dig-2); 1611 1612 assert(feqrel(real.min_normal/8,real.min_normal/17)==3); 1613 1614 // Numbers that are close 1615 assert(feqrel(0x1.Bp+84, 0x1.B8p+84)==5); 1616 assert(feqrel(0x1.8p+10, 0x1.Cp+10)==2); 1617 assert(feqrel(1.5*(1-real.epsilon), 1.0L)==2); 1618 assert(feqrel(1.5, 1.0)==1); 1619 assert(feqrel(2*(1-real.epsilon), 1.0L)==1); 1620 1621 // Factors of 2 1622 assert(feqrel(real.max,real.infinity)==0); 1623 assert(feqrel(2*(1-real.epsilon), 1.0L)==1); 1624 assert(feqrel(1.0, 2.0)==0); 1625 assert(feqrel(4.0, 1.0)==0); 1626 1627 // Extreme inequality 1628 assert(feqrel(real.nan,real.nan)==0); 1629 assert(feqrel(0.0L,-real.nan)==0); 1630 assert(feqrel(real.nan,real.infinity)==0); 1631 assert(feqrel(real.infinity,-real.infinity)==0); 1632 assert(feqrel(-real.max,real.infinity)==0); 1633 assert(feqrel(real.max,-real.max)==0); 1634 1635 // floats 1636 assert(feqrel(2.1f, 2.1f)==float.mant_dig); 1637 assert(feqrel(1.5f, 1.0f)==1); 1638 } 1639 } 1640 1641 /********************************* 1642 * Return 1 if sign bit of e is set, 0 if not. 1643 */ 1644 1645 int signbit(real x) 1646 { 1647 return ((cast(ubyte *)&x)[floatTraits!(real).SIGNPOS_BYTE] & 0x80) != 0; 1648 } 1649 1650 debug(UnitTest) { 1651 unittest 1652 { 1653 assert(!signbit(float.nan)); 1654 assert(signbit(-float.nan)); 1655 assert(!signbit(168.1234)); 1656 assert(signbit(-168.1234)); 1657 assert(!signbit(0.0)); 1658 assert(signbit(-0.0)); 1659 } 1660 } 1661 1662 1663 /********************************* 1664 * Return a value composed of to with from's sign bit. 1665 */ 1666 1667 real copysign(real to, real from) 1668 { 1669 ubyte* pto = cast(ubyte *)&to; 1670 ubyte* pfrom = cast(ubyte *)&from; 1671 1672 alias floatTraits!(real) F; 1673 pto[F.SIGNPOS_BYTE] &= 0x7F; 1674 pto[F.SIGNPOS_BYTE] |= pfrom[F.SIGNPOS_BYTE] & 0x80; 1675 return to; 1676 } 1677 1678 debug(UnitTest) { 1679 unittest 1680 { 1681 real e; 1682 1683 e = copysign(21, 23.8); 1684 assert(e == 21); 1685 1686 e = copysign(-21, 23.8); 1687 assert(e == 21); 1688 1689 e = copysign(21, -23.8); 1690 assert(e == -21); 1691 1692 e = copysign(-21, -23.8); 1693 assert(e == -21); 1694 1695 e = copysign(real.nan, -23.8); 1696 assert(isNaN(e) && signbit(e)); 1697 } 1698 } 1699 1700 /** Return the value that lies halfway between x and y on the IEEE number line. 1701 * 1702 * Formally, the result is the arithmetic mean of the binary significands of x 1703 * and y, multiplied by the geometric mean of the binary exponents of x and y. 1704 * x and y must have the same sign, and must not be NaN. 1705 * Note: this function is useful for ensuring O(log n) behaviour in algorithms 1706 * involving a 'binary chop'. 1707 * 1708 * Special cases: 1709 * If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value 1710 * is the arithmetic mean (x + y) / 2. 1711 * If x and y are even powers of 2, the return value is the geometric mean, 1712 * ieeeMean(x, y) = sqrt(x * y). 1713 * 1714 */ 1715 T ieeeMean(T)(T x, T y) 1716 in { 1717 // both x and y must have the same sign, and must not be NaN. 1718 assert(signbit(x) == signbit(y)); 1719 assert(!isNaN(x) && !isNaN(y)); 1720 } 1721 body { 1722 // Runtime behaviour for contract violation: 1723 // If signs are opposite, or one is a NaN, return 0. 1724 if (!((x>=0 && y>=0) || (x<=0 && y<=0))) return 0.0; 1725 1726 // The implementation is simple: cast x and y to integers, 1727 // average them (avoiding overflow), and cast the result back to a floating-point number. 1728 1729 alias floatTraits!(real) F; 1730 T u; 1731 static if (T.mant_dig==64) { // real80 1732 // There's slight additional complexity because they are actually 1733 // 79-bit reals... 1734 ushort *ue = cast(ushort *)&u; 1735 ulong *ul = cast(ulong *)&u; 1736 ushort *xe = cast(ushort *)&x; 1737 ulong *xl = cast(ulong *)&x; 1738 ushort *ye = cast(ushort *)&y; 1739 ulong *yl = cast(ulong *)&y; 1740 // Ignore the useless implicit bit. (Bonus: this prevents overflows) 1741 ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL); 1742 1743 ushort e = cast(ushort)((xe[F.EXPPOS_SHORT] & 0x7FFF) + (ye[F.EXPPOS_SHORT] & 0x7FFF)); 1744 if (m & 0x8000_0000_0000_0000L) { 1745 ++e; 1746 m &= 0x7FFF_FFFF_FFFF_FFFFL; 1747 } 1748 // Now do a multi-byte right shift 1749 uint c = e & 1; // carry 1750 e >>= 1; 1751 m >>>= 1; 1752 if (c) m |= 0x4000_0000_0000_0000L; // shift carry into significand 1753 if (e) *ul = m | 0x8000_0000_0000_0000L; // set implicit bit... 1754 else *ul = m; // ... unless exponent is 0 (denormal or zero). 1755 ue[4]= e | (xe[F.EXPPOS_SHORT]& F.SIGNMASK); // restore sign bit 1756 } else static if(T.mant_dig == 113) { //quadruple 1757 // This would be trivial if 'ucent' were implemented... 1758 ulong *ul = cast(ulong *)&u; 1759 ulong *xl = cast(ulong *)&x; 1760 ulong *yl = cast(ulong *)&y; 1761 // Multi-byte add, then multi-byte right shift. 1762 ulong mh = ((xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL) 1763 + (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL)); 1764 // Discard the lowest bit (to avoid overflow) 1765 ulong ml = (xl[MANTISSA_LSB]>>>1) + (yl[MANTISSA_LSB]>>>1); 1766 // add the lowest bit back in, if necessary. 1767 if (xl[MANTISSA_LSB] & yl[MANTISSA_LSB] & 1) { 1768 ++ml; 1769 if (ml==0) ++mh; 1770 } 1771 mh >>>=1; 1772 ul[MANTISSA_MSB] = mh | (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000); 1773 ul[MANTISSA_LSB] = ml; 1774 } else static if (T.mant_dig == double.mant_dig) { 1775 ulong *ul = cast(ulong *)&u; 1776 ulong *xl = cast(ulong *)&x; 1777 ulong *yl = cast(ulong *)&y; 1778 ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1; 1779 m |= ((*xl) & 0x8000_0000_0000_0000L); 1780 *ul = m; 1781 } else static if (T.mant_dig == float.mant_dig) { 1782 uint *ul = cast(uint *)&u; 1783 uint *xl = cast(uint *)&x; 1784 uint *yl = cast(uint *)&y; 1785 uint m = (((*xl) & 0x7FFF_FFFF) + ((*yl) & 0x7FFF_FFFF)) >>> 1; 1786 m |= ((*xl) & 0x8000_0000); 1787 *ul = m; 1788 } else { 1789 assert(0, "Not implemented"); 1790 } 1791 return u; 1792 } 1793 1794 debug(UnitTest) { 1795 unittest { 1796 assert(ieeeMean(-0.0,-1e-20)<0); 1797 assert(ieeeMean(0.0,1e-20)>0); 1798 1799 assert(ieeeMean(1.0L,4.0L)==2L); 1800 assert(ieeeMean(2.0*1.013,8.0*1.013)==4*1.013); 1801 assert(ieeeMean(-1.0L,-4.0L)==-2L); 1802 assert(ieeeMean(-1.0,-4.0)==-2); 1803 assert(ieeeMean(-1.0f,-4.0f)==-2f); 1804 assert(ieeeMean(-1.0,-2.0)==-1.5); 1805 assert(ieeeMean(-1*(1+8*real.epsilon),-2*(1+8*real.epsilon))==-1.5*(1+5*real.epsilon)); 1806 assert(ieeeMean(0x1p60,0x1p-10)==0x1p25); 1807 static if (real.mant_dig==64) { // x87, 80-bit reals 1808 assert(ieeeMean(1.0L,real.infinity)==0x1p8192L); 1809 assert(ieeeMean(0.0L,real.infinity)==1.5); 1810 } 1811 assert(ieeeMean(0.5*real.min_normal*(1-4*real.epsilon),0.5*real.min_normal)==0.5*real.min_normal*(1-2*real.epsilon)); 1812 } 1813 } 1814 1815 // Functions for NaN payloads 1816 /* 1817 * A 'payload' can be stored in the significand of a $(NAN). One bit is required 1818 * to distinguish between a quiet and a signalling $(NAN). This leaves 22 bits 1819 * of payload for a float; 51 bits for a double; 62 bits for an 80-bit real; 1820 * and 111 bits for a 128-bit quad. 1821 */ 1822 /** 1823 * Create a $(NAN), storing an integer inside the payload. 1824 * 1825 * For 80-bit or 128-bit reals, the largest possible payload is 0x3FFF_FFFF_FFFF_FFFF. 1826 * For doubles, it is 0x3_FFFF_FFFF_FFFF. 1827 * For floats, it is 0x3F_FFFF. 1828 */ 1829 real NaN(ulong payload) 1830 { 1831 static if (real.mant_dig == 64) { //real80 1832 ulong v = 3; // implied bit = 1, quiet bit = 1 1833 } else { 1834 ulong v = 2; // no implied bit. quiet bit = 1 1835 } 1836 1837 ulong a = payload; 1838 1839 // 22 Float bits 1840 ulong w = a & 0x3F_FFFF; 1841 a -= w; 1842 1843 v <<=22; 1844 v |= w; 1845 a >>=22; 1846 1847 // 29 Double bits 1848 v <<=29; 1849 w = a & 0xFFF_FFFF; 1850 v |= w; 1851 a -= w; 1852 a >>=29; 1853 1854 static if (real.mant_dig == 53) { // double 1855 v |=0x7FF0_0000_0000_0000; 1856 real x; 1857 * cast(ulong *)(&x) = v; 1858 return x; 1859 } else { 1860 v <<=11; 1861 a &= 0x7FF; 1862 v |= a; 1863 real x = real.nan; 1864 // Extended real bits 1865 static if (real.mant_dig==113) { //quadruple 1866 v<<=1; // there's no implicit bit 1867 version(LittleEndian) { 1868 *cast(ulong*)(6+cast(ubyte*)(&x)) = v; 1869 } else { 1870 *cast(ulong*)(2+cast(ubyte*)(&x)) = v; 1871 } 1872 } else { // real80 1873 * cast(ulong *)(&x) = v; 1874 } 1875 return x; 1876 } 1877 } 1878 1879 /** 1880 * Extract an integral payload from a $(NAN). 1881 * 1882 * Returns: 1883 * the integer payload as a ulong. 1884 * 1885 * For 80-bit or 128-bit reals, the largest possible payload is 0x3FFF_FFFF_FFFF_FFFF. 1886 * For doubles, it is 0x3_FFFF_FFFF_FFFF. 1887 * For floats, it is 0x3F_FFFF. 1888 */ 1889 ulong getNaNPayload(real x) 1890 { 1891 assert(isNaN(x)); 1892 static if (real.mant_dig == 53) { 1893 ulong m = *cast(ulong *)(&x); 1894 // Make it look like an 80-bit significand. 1895 // Skip exponent, and quiet bit 1896 m &= 0x0007_FFFF_FFFF_FFFF; 1897 m <<= 10; 1898 } else static if (real.mant_dig==113) { // quadruple 1899 version(LittleEndian) { 1900 ulong m = *cast(ulong*)(6+cast(ubyte*)(&x)); 1901 } else { 1902 ulong m = *cast(ulong*)(2+cast(ubyte*)(&x)); 1903 } 1904 m>>=1; // there's no implicit bit 1905 } else { 1906 ulong m = *cast(ulong *)(&x); 1907 } 1908 // ignore implicit bit and quiet bit 1909 ulong f = m & 0x3FFF_FF00_0000_0000L; 1910 ulong w = f >>> 40; 1911 w |= (m & 0x00FF_FFFF_F800L) << (22 - 11); 1912 w |= (m & 0x7FF) << 51; 1913 return w; 1914 } 1915 1916 debug(UnitTest) { 1917 unittest { 1918 real nan4 = NaN(0x789_ABCD_EF12_3456); 1919 static if (real.mant_dig == 64 || real.mant_dig==113) { 1920 assert (getNaNPayload(nan4) == 0x789_ABCD_EF12_3456); 1921 } else { 1922 assert (getNaNPayload(nan4) == 0x1_ABCD_EF12_3456); 1923 } 1924 double nan5 = nan4; 1925 assert (getNaNPayload(nan5) == 0x1_ABCD_EF12_3456); 1926 float nan6 = nan4; 1927 assert (getNaNPayload(nan6) == 0x12_3456); 1928 nan4 = NaN(0xFABCD); 1929 assert (getNaNPayload(nan4) == 0xFABCD); 1930 nan6 = nan4; 1931 assert (getNaNPayload(nan6) == 0xFABCD); 1932 nan5 = NaN(0x100_0000_0000_3456); 1933 assert(getNaNPayload(nan5) == 0x0000_0000_3456); 1934 } 1935 } 1936