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