1 // Written in the D programming language. 2 3 /** 4 This is a submodule of $(MREF std, math). 5 6 It contains several functions for rounding floating point numbers. 7 8 Copyright: Copyright The D Language Foundation 2000 - 2011. 9 D implementations of floor, ceil, and lrint functions are based on the 10 CEPHES math library, which is Copyright (C) 2001 Stephen L. Moshier 11 $(LT)steve@moshier.net$(GT) and are incorporated herein by permission 12 of the author. The author reserves the right to distribute this 13 material elsewhere under different copying permissions. 14 These modifications are distributed here under the following terms: 15 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). 16 Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston, 17 Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger 18 Source: $(PHOBOSSRC std/math/rounding.d) 19 */ 20 21 module std.math.rounding; 22 23 static import core.math; 24 static import core.stdc.math; 25 26 import std.traits : isFloatingPoint, isIntegral, Unqual; 27 28 version (D_InlineAsm_X86) version = InlineAsm_X86_Any; 29 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any; 30 31 version (InlineAsm_X86_Any) version = InlineAsm_X87; 32 version (InlineAsm_X87) 33 { 34 static assert(real.mant_dig == 64); 35 version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC; 36 } 37 38 /************************************** 39 * Returns the value of x rounded upward to the next integer 40 * (toward positive infinity). 41 */ 42 real ceil(real x) @trusted pure nothrow @nogc 43 { 44 version (InlineAsm_X87_MSVC) 45 { 46 version (X86_64) 47 { 48 asm pure nothrow @nogc 49 { 50 naked ; 51 fld real ptr [RCX] ; 52 fstcw 8[RSP] ; 53 mov AL,9[RSP] ; 54 mov DL,AL ; 55 and AL,0xC3 ; 56 or AL,0x08 ; // round to +infinity 57 mov 9[RSP],AL ; 58 fldcw 8[RSP] ; 59 frndint ; 60 mov 9[RSP],DL ; 61 fldcw 8[RSP] ; 62 ret ; 63 } 64 } 65 else 66 { 67 short cw; 68 asm pure nothrow @nogc 69 { 70 fld x ; 71 fstcw cw ; 72 mov AL,byte ptr cw+1 ; 73 mov DL,AL ; 74 and AL,0xC3 ; 75 or AL,0x08 ; // round to +infinity 76 mov byte ptr cw+1,AL ; 77 fldcw cw ; 78 frndint ; 79 mov byte ptr cw+1,DL ; 80 fldcw cw ; 81 } 82 } 83 } 84 else 85 { 86 import std.math.traits : isInfinity, isNaN; 87 88 // Special cases. 89 if (isNaN(x) || isInfinity(x)) 90 return x; 91 92 real y = floorImpl(x); 93 if (y < x) 94 y += 1.0; 95 96 return y; 97 } 98 } 99 100 /// 101 @safe pure nothrow @nogc unittest 102 { 103 import std.math.traits : isNaN; 104 105 assert(ceil(+123.456L) == +124); 106 assert(ceil(-123.456L) == -123); 107 assert(ceil(-1.234L) == -1); 108 assert(ceil(-0.123L) == 0); 109 assert(ceil(0.0L) == 0); 110 assert(ceil(+0.123L) == 1); 111 assert(ceil(+1.234L) == 2); 112 assert(ceil(real.infinity) == real.infinity); 113 assert(isNaN(ceil(real.nan))); 114 assert(isNaN(ceil(real.init))); 115 } 116 117 /// ditto 118 double ceil(double x) @trusted pure nothrow @nogc 119 { 120 import std.math.traits : isInfinity, isNaN; 121 122 // Special cases. 123 if (isNaN(x) || isInfinity(x)) 124 return x; 125 126 double y = floorImpl(x); 127 if (y < x) 128 y += 1.0; 129 130 return y; 131 } 132 133 @safe pure nothrow @nogc unittest 134 { 135 import std.math.traits : isNaN; 136 137 assert(ceil(+123.456) == +124); 138 assert(ceil(-123.456) == -123); 139 assert(ceil(-1.234) == -1); 140 assert(ceil(-0.123) == 0); 141 assert(ceil(0.0) == 0); 142 assert(ceil(+0.123) == 1); 143 assert(ceil(+1.234) == 2); 144 assert(ceil(double.infinity) == double.infinity); 145 assert(isNaN(ceil(double.nan))); 146 assert(isNaN(ceil(double.init))); 147 } 148 149 /// ditto 150 float ceil(float x) @trusted pure nothrow @nogc 151 { 152 import std.math.traits : isInfinity, isNaN; 153 154 // Special cases. 155 if (isNaN(x) || isInfinity(x)) 156 return x; 157 158 float y = floorImpl(x); 159 if (y < x) 160 y += 1.0; 161 162 return y; 163 } 164 165 @safe pure nothrow @nogc unittest 166 { 167 import std.math.traits : isNaN; 168 169 assert(ceil(+123.456f) == +124); 170 assert(ceil(-123.456f) == -123); 171 assert(ceil(-1.234f) == -1); 172 assert(ceil(-0.123f) == 0); 173 assert(ceil(0.0f) == 0); 174 assert(ceil(+0.123f) == 1); 175 assert(ceil(+1.234f) == 2); 176 assert(ceil(float.infinity) == float.infinity); 177 assert(isNaN(ceil(float.nan))); 178 assert(isNaN(ceil(float.init))); 179 } 180 181 /************************************** 182 * Returns the value of x rounded downward to the next integer 183 * (toward negative infinity). 184 */ 185 real floor(real x) @trusted pure nothrow @nogc 186 { 187 version (InlineAsm_X87_MSVC) 188 { 189 version (X86_64) 190 { 191 asm pure nothrow @nogc 192 { 193 naked ; 194 fld real ptr [RCX] ; 195 fstcw 8[RSP] ; 196 mov AL,9[RSP] ; 197 mov DL,AL ; 198 and AL,0xC3 ; 199 or AL,0x04 ; // round to -infinity 200 mov 9[RSP],AL ; 201 fldcw 8[RSP] ; 202 frndint ; 203 mov 9[RSP],DL ; 204 fldcw 8[RSP] ; 205 ret ; 206 } 207 } 208 else 209 { 210 short cw; 211 asm pure nothrow @nogc 212 { 213 fld x ; 214 fstcw cw ; 215 mov AL,byte ptr cw+1 ; 216 mov DL,AL ; 217 and AL,0xC3 ; 218 or AL,0x04 ; // round to -infinity 219 mov byte ptr cw+1,AL ; 220 fldcw cw ; 221 frndint ; 222 mov byte ptr cw+1,DL ; 223 fldcw cw ; 224 } 225 } 226 } 227 else 228 { 229 import std.math.traits : isInfinity, isNaN; 230 231 // Special cases. 232 if (isNaN(x) || isInfinity(x) || x == 0.0) 233 return x; 234 235 return floorImpl(x); 236 } 237 } 238 239 /// 240 @safe pure nothrow @nogc unittest 241 { 242 import std.math.traits : isNaN; 243 244 assert(floor(+123.456L) == +123); 245 assert(floor(-123.456L) == -124); 246 assert(floor(+123.0L) == +123); 247 assert(floor(-124.0L) == -124); 248 assert(floor(-1.234L) == -2); 249 assert(floor(-0.123L) == -1); 250 assert(floor(0.0L) == 0); 251 assert(floor(+0.123L) == 0); 252 assert(floor(+1.234L) == 1); 253 assert(floor(real.infinity) == real.infinity); 254 assert(isNaN(floor(real.nan))); 255 assert(isNaN(floor(real.init))); 256 } 257 258 /// ditto 259 double floor(double x) @trusted pure nothrow @nogc 260 { 261 import std.math.traits : isInfinity, isNaN; 262 263 // Special cases. 264 if (isNaN(x) || isInfinity(x) || x == 0.0) 265 return x; 266 267 return floorImpl(x); 268 } 269 270 @safe pure nothrow @nogc unittest 271 { 272 import std.math.traits : isNaN; 273 274 assert(floor(+123.456) == +123); 275 assert(floor(-123.456) == -124); 276 assert(floor(+123.0) == +123); 277 assert(floor(-124.0) == -124); 278 assert(floor(-1.234) == -2); 279 assert(floor(-0.123) == -1); 280 assert(floor(0.0) == 0); 281 assert(floor(+0.123) == 0); 282 assert(floor(+1.234) == 1); 283 assert(floor(double.infinity) == double.infinity); 284 assert(isNaN(floor(double.nan))); 285 assert(isNaN(floor(double.init))); 286 } 287 288 /// ditto 289 float floor(float x) @trusted pure nothrow @nogc 290 { 291 import std.math.traits : isInfinity, isNaN; 292 293 // Special cases. 294 if (isNaN(x) || isInfinity(x) || x == 0.0) 295 return x; 296 297 return floorImpl(x); 298 } 299 300 @safe pure nothrow @nogc unittest 301 { 302 import std.math.traits : isNaN; 303 304 assert(floor(+123.456f) == +123); 305 assert(floor(-123.456f) == -124); 306 assert(floor(+123.0f) == +123); 307 assert(floor(-124.0f) == -124); 308 assert(floor(-1.234f) == -2); 309 assert(floor(-0.123f) == -1); 310 assert(floor(0.0f) == 0); 311 assert(floor(+0.123f) == 0); 312 assert(floor(+1.234f) == 1); 313 assert(floor(float.infinity) == float.infinity); 314 assert(isNaN(floor(float.nan))); 315 assert(isNaN(floor(float.init))); 316 } 317 318 // https://issues.dlang.org/show_bug.cgi?id=6381 319 // floor/ceil should be usable in pure function. 320 @safe pure nothrow unittest 321 { 322 auto x = floor(1.2); 323 auto y = ceil(1.2); 324 } 325 326 /** 327 * Round `val` to a multiple of `unit`. `rfunc` specifies the rounding 328 * function to use; by default this is `rint`, which uses the current 329 * rounding mode. 330 */ 331 Unqual!F quantize(alias rfunc = rint, F)(const F val, const F unit) 332 if (is(typeof(rfunc(F.init)) : F) && isFloatingPoint!F) 333 { 334 import std.math.traits : isInfinity; 335 336 typeof(return) ret = val; 337 if (unit != 0) 338 { 339 const scaled = val / unit; 340 if (!scaled.isInfinity) 341 ret = rfunc(scaled) * unit; 342 } 343 return ret; 344 } 345 346 /// 347 @safe pure nothrow @nogc unittest 348 { 349 import std.math.operations : isClose; 350 351 assert(isClose(12345.6789L.quantize(0.01L), 12345.68L)); 352 assert(isClose(12345.6789L.quantize!floor(0.01L), 12345.67L)); 353 assert(isClose(12345.6789L.quantize(22.0L), 12342.0L)); 354 } 355 356 /// 357 @safe pure nothrow @nogc unittest 358 { 359 import std.math.operations : isClose; 360 import std.math.traits : isNaN; 361 362 assert(isClose(12345.6789L.quantize(0), 12345.6789L)); 363 assert(12345.6789L.quantize(real.infinity).isNaN); 364 assert(12345.6789L.quantize(real.nan).isNaN); 365 assert(real.infinity.quantize(0.01L) == real.infinity); 366 assert(real.infinity.quantize(real.nan).isNaN); 367 assert(real.nan.quantize(0.01L).isNaN); 368 assert(real.nan.quantize(real.infinity).isNaN); 369 assert(real.nan.quantize(real.nan).isNaN); 370 } 371 372 /** 373 * Round `val` to a multiple of `pow(base, exp)`. `rfunc` specifies the 374 * rounding function to use; by default this is `rint`, which uses the 375 * current rounding mode. 376 */ 377 Unqual!F quantize(real base, alias rfunc = rint, F, E)(const F val, const E exp) 378 if (is(typeof(rfunc(F.init)) : F) && isFloatingPoint!F && isIntegral!E) 379 { 380 import std.math.exponential : pow; 381 382 // TODO: Compile-time optimization for power-of-two bases? 383 return quantize!rfunc(val, pow(cast(F) base, exp)); 384 } 385 386 /// ditto 387 Unqual!F quantize(real base, long exp = 1, alias rfunc = rint, F)(const F val) 388 if (is(typeof(rfunc(F.init)) : F) && isFloatingPoint!F) 389 { 390 import std.math.exponential : pow; 391 392 enum unit = cast(F) pow(base, exp); 393 return quantize!rfunc(val, unit); 394 } 395 396 /// 397 @safe pure nothrow @nogc unittest 398 { 399 import std.math.operations : isClose; 400 401 assert(isClose(12345.6789L.quantize!10(-2), 12345.68L)); 402 assert(isClose(12345.6789L.quantize!(10, -2), 12345.68L)); 403 assert(isClose(12345.6789L.quantize!(10, floor)(-2), 12345.67L)); 404 assert(isClose(12345.6789L.quantize!(10, -2, floor), 12345.67L)); 405 406 assert(isClose(12345.6789L.quantize!22(1), 12342.0L)); 407 assert(isClose(12345.6789L.quantize!22, 12342.0L)); 408 } 409 410 @safe pure nothrow @nogc unittest 411 { 412 import std.math.exponential : log10, pow; 413 import std.math.operations : isClose; 414 import std.meta : AliasSeq; 415 416 static foreach (F; AliasSeq!(real, double, float)) 417 {{ 418 const maxL10 = cast(int) F.max.log10.floor; 419 const maxR10 = pow(cast(F) 10, maxL10); 420 assert(isClose((cast(F) 0.9L * maxR10).quantize!10(maxL10), maxR10)); 421 assert(isClose((cast(F)-0.9L * maxR10).quantize!10(maxL10), -maxR10)); 422 423 assert(F.max.quantize(F.min_normal) == F.max); 424 assert((-F.max).quantize(F.min_normal) == -F.max); 425 assert(F.min_normal.quantize(F.max) == 0); 426 assert((-F.min_normal).quantize(F.max) == 0); 427 assert(F.min_normal.quantize(F.min_normal) == F.min_normal); 428 assert((-F.min_normal).quantize(F.min_normal) == -F.min_normal); 429 }} 430 } 431 432 /****************************************** 433 * Rounds x to the nearest integer value, using the current rounding 434 * mode. 435 * 436 * Unlike the rint functions, nearbyint does not raise the 437 * FE_INEXACT exception. 438 */ 439 pragma(inline, true) 440 real nearbyint(real x) @safe pure nothrow @nogc 441 { 442 return core.stdc.math.nearbyintl(x); 443 } 444 445 /// 446 @safe pure unittest 447 { 448 import std.math.traits : isNaN; 449 450 assert(nearbyint(0.4) == 0); 451 assert(nearbyint(0.5) == 0); 452 assert(nearbyint(0.6) == 1); 453 assert(nearbyint(100.0) == 100); 454 455 assert(isNaN(nearbyint(real.nan))); 456 assert(nearbyint(real.infinity) == real.infinity); 457 assert(nearbyint(-real.infinity) == -real.infinity); 458 } 459 460 /********************************** 461 * Rounds x to the nearest integer value, using the current rounding 462 * mode. 463 * 464 * If the return value is not equal to x, the FE_INEXACT 465 * exception is raised. 466 * 467 * $(LREF nearbyint) performs the same operation, but does 468 * not set the FE_INEXACT exception. 469 */ 470 pragma(inline, true) 471 real rint(real x) @safe pure nothrow @nogc 472 { 473 return core.math.rint(x); 474 } 475 ///ditto 476 pragma(inline, true) 477 double rint(double x) @safe pure nothrow @nogc 478 { 479 return core.math.rint(x); 480 } 481 ///ditto 482 pragma(inline, true) 483 float rint(float x) @safe pure nothrow @nogc 484 { 485 return core.math.rint(x); 486 } 487 488 /// 489 @safe unittest 490 { 491 import std.math.traits : isNaN; 492 493 version (IeeeFlagsSupport) resetIeeeFlags(); 494 assert(rint(0.4) == 0); 495 version (IeeeFlagsSupport) assert(ieeeFlags.inexact); 496 497 assert(rint(0.5) == 0); 498 assert(rint(0.6) == 1); 499 assert(rint(100.0) == 100); 500 501 assert(isNaN(rint(real.nan))); 502 assert(rint(real.infinity) == real.infinity); 503 assert(rint(-real.infinity) == -real.infinity); 504 } 505 506 @safe unittest 507 { 508 real function(real) print = &rint; 509 assert(print != null); 510 } 511 512 /*************************************** 513 * Rounds x to the nearest integer value, using the current rounding 514 * mode. 515 * 516 * This is generally the fastest method to convert a floating-point number 517 * to an integer. Note that the results from this function 518 * depend on the rounding mode, if the fractional part of x is exactly 0.5. 519 * If using the default rounding mode (ties round to even integers) 520 * lrint(4.5) == 4, lrint(5.5)==6. 521 */ 522 long lrint(real x) @trusted pure nothrow @nogc 523 { 524 version (InlineAsm_X87) 525 { 526 version (Win64) 527 { 528 asm pure nothrow @nogc 529 { 530 naked; 531 fld real ptr [RCX]; 532 fistp qword ptr 8[RSP]; 533 mov RAX,8[RSP]; 534 ret; 535 } 536 } 537 else 538 { 539 long n; 540 asm pure nothrow @nogc 541 { 542 fld x; 543 fistp n; 544 } 545 return n; 546 } 547 } 548 else 549 { 550 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB; 551 552 alias F = floatTraits!(real); 553 static if (F.realFormat == RealFormat.ieeeDouble) 554 { 555 long result; 556 557 // Rounding limit when casting from real(double) to ulong. 558 enum real OF = 4.50359962737049600000E15L; 559 560 uint* vi = cast(uint*)(&x); 561 562 // Find the exponent and sign 563 uint msb = vi[MANTISSA_MSB]; 564 uint lsb = vi[MANTISSA_LSB]; 565 int exp = ((msb >> 20) & 0x7ff) - 0x3ff; 566 const int sign = msb >> 31; 567 msb &= 0xfffff; 568 msb |= 0x100000; 569 570 if (exp < 63) 571 { 572 if (exp >= 52) 573 result = (cast(long) msb << (exp - 20)) | (lsb << (exp - 52)); 574 else 575 { 576 // Adjust x and check result. 577 const real j = sign ? -OF : OF; 578 x = (j + x) - j; 579 msb = vi[MANTISSA_MSB]; 580 lsb = vi[MANTISSA_LSB]; 581 exp = ((msb >> 20) & 0x7ff) - 0x3ff; 582 msb &= 0xfffff; 583 msb |= 0x100000; 584 585 if (exp < 0) 586 result = 0; 587 else if (exp < 20) 588 result = cast(long) msb >> (20 - exp); 589 else if (exp == 20) 590 result = cast(long) msb; 591 else 592 result = (cast(long) msb << (exp - 20)) | (lsb >> (52 - exp)); 593 } 594 } 595 else 596 { 597 // It is left implementation defined when the number is too large. 598 return cast(long) x; 599 } 600 601 return sign ? -result : result; 602 } 603 else static if (F.realFormat == RealFormat.ieeeExtended || 604 F.realFormat == RealFormat.ieeeExtended53) 605 { 606 long result; 607 608 // Rounding limit when casting from real(80-bit) to ulong. 609 static if (F.realFormat == RealFormat.ieeeExtended) 610 enum real OF = 9.22337203685477580800E18L; 611 else 612 enum real OF = 4.50359962737049600000E15L; 613 614 ushort* vu = cast(ushort*)(&x); 615 uint* vi = cast(uint*)(&x); 616 617 // Find the exponent and sign 618 int exp = (vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff; 619 const int sign = (vu[F.EXPPOS_SHORT] >> 15) & 1; 620 621 if (exp < 63) 622 { 623 // Adjust x and check result. 624 const real j = sign ? -OF : OF; 625 x = (j + x) - j; 626 exp = (vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff; 627 628 version (LittleEndian) 629 { 630 if (exp < 0) 631 result = 0; 632 else if (exp <= 31) 633 result = vi[1] >> (31 - exp); 634 else 635 result = (cast(long) vi[1] << (exp - 31)) | (vi[0] >> (63 - exp)); 636 } 637 else 638 { 639 if (exp < 0) 640 result = 0; 641 else if (exp <= 31) 642 result = vi[1] >> (31 - exp); 643 else 644 result = (cast(long) vi[1] << (exp - 31)) | (vi[2] >> (63 - exp)); 645 } 646 } 647 else 648 { 649 // It is left implementation defined when the number is too large 650 // to fit in a 64bit long. 651 return cast(long) x; 652 } 653 654 return sign ? -result : result; 655 } 656 else static if (F.realFormat == RealFormat.ieeeQuadruple) 657 { 658 const vu = cast(ushort*)(&x); 659 660 // Find the exponent and sign 661 const sign = (vu[F.EXPPOS_SHORT] >> 15) & 1; 662 if ((vu[F.EXPPOS_SHORT] & F.EXPMASK) - (F.EXPBIAS + 1) > 63) 663 { 664 // The result is left implementation defined when the number is 665 // too large to fit in a 64 bit long. 666 return cast(long) x; 667 } 668 669 // Force rounding of lower bits according to current rounding 670 // mode by adding ±2^-112 and subtracting it again. 671 enum OF = 5.19229685853482762853049632922009600E33L; 672 const j = sign ? -OF : OF; 673 x = (j + x) - j; 674 675 const exp = (vu[F.EXPPOS_SHORT] & F.EXPMASK) - (F.EXPBIAS + 1); 676 const implicitOne = 1UL << 48; 677 auto vl = cast(ulong*)(&x); 678 vl[MANTISSA_MSB] &= implicitOne - 1; 679 vl[MANTISSA_MSB] |= implicitOne; 680 681 long result; 682 683 if (exp < 0) 684 result = 0; 685 else if (exp <= 48) 686 result = vl[MANTISSA_MSB] >> (48 - exp); 687 else 688 result = (vl[MANTISSA_MSB] << (exp - 48)) | (vl[MANTISSA_LSB] >> (112 - exp)); 689 690 return sign ? -result : result; 691 } 692 else 693 { 694 static assert(false, "real type not supported by lrint()"); 695 } 696 } 697 } 698 699 /// 700 @safe pure nothrow @nogc unittest 701 { 702 assert(lrint(4.5) == 4); 703 assert(lrint(5.5) == 6); 704 assert(lrint(-4.5) == -4); 705 assert(lrint(-5.5) == -6); 706 707 assert(lrint(int.max - 0.5) == 2147483646L); 708 assert(lrint(int.max + 0.5) == 2147483648L); 709 assert(lrint(int.min - 0.5) == -2147483648L); 710 assert(lrint(int.min + 0.5) == -2147483648L); 711 } 712 713 static if (real.mant_dig >= long.sizeof * 8) 714 { 715 @safe pure nothrow @nogc unittest 716 { 717 assert(lrint(long.max - 1.5L) == long.max - 1); 718 assert(lrint(long.max - 0.5L) == long.max - 1); 719 assert(lrint(long.min + 0.5L) == long.min); 720 assert(lrint(long.min + 1.5L) == long.min + 2); 721 } 722 } 723 724 /******************************************* 725 * Return the value of x rounded to the nearest integer. 726 * If the fractional part of x is exactly 0.5, the return value is 727 * rounded away from zero. 728 * 729 * Returns: 730 * A `real`. 731 */ 732 auto round(real x) @trusted nothrow @nogc 733 { 734 version (CRuntime_Microsoft) 735 { 736 import std.math.hardware : FloatingPointControl; 737 738 auto old = FloatingPointControl.getControlState(); 739 FloatingPointControl.setControlState( 740 (old & (-1 - FloatingPointControl.roundingMask)) | FloatingPointControl.roundToZero 741 ); 742 x = core.math.rint((x >= 0) ? x + 0.5 : x - 0.5); 743 FloatingPointControl.setControlState(old); 744 return x; 745 } 746 else 747 { 748 return core.stdc.math.roundl(x); 749 } 750 } 751 752 /// 753 @safe nothrow @nogc unittest 754 { 755 assert(round(4.5) == 5); 756 assert(round(5.4) == 5); 757 assert(round(-4.5) == -5); 758 assert(round(-5.1) == -5); 759 } 760 761 // assure purity on Posix 762 version (Posix) 763 { 764 @safe pure nothrow @nogc unittest 765 { 766 assert(round(4.5) == 5); 767 } 768 } 769 770 /********************************************** 771 * Return the value of x rounded to the nearest integer. 772 * 773 * If the fractional part of x is exactly 0.5, the return value is rounded 774 * away from zero. 775 * 776 * $(BLUE This function is not implemented for Digital Mars C runtime.) 777 */ 778 long lround(real x) @trusted nothrow @nogc 779 { 780 version (CRuntime_DigitalMars) 781 assert(0, "lround not implemented"); 782 else 783 return core.stdc.math.llroundl(x); 784 } 785 786 /// 787 @safe nothrow @nogc unittest 788 { 789 version (CRuntime_DigitalMars) {} 790 else 791 { 792 assert(lround(0.49) == 0); 793 assert(lround(0.5) == 1); 794 assert(lround(1.5) == 2); 795 } 796 } 797 798 /** 799 Returns the integer portion of x, dropping the fractional portion. 800 This is also known as "chop" rounding. 801 `pure` on all platforms. 802 */ 803 real trunc(real x) @trusted nothrow @nogc pure 804 { 805 version (InlineAsm_X87_MSVC) 806 { 807 version (X86_64) 808 { 809 asm pure nothrow @nogc 810 { 811 naked ; 812 fld real ptr [RCX] ; 813 fstcw 8[RSP] ; 814 mov AL,9[RSP] ; 815 mov DL,AL ; 816 and AL,0xC3 ; 817 or AL,0x0C ; // round to 0 818 mov 9[RSP],AL ; 819 fldcw 8[RSP] ; 820 frndint ; 821 mov 9[RSP],DL ; 822 fldcw 8[RSP] ; 823 ret ; 824 } 825 } 826 else 827 { 828 short cw; 829 asm pure nothrow @nogc 830 { 831 fld x ; 832 fstcw cw ; 833 mov AL,byte ptr cw+1 ; 834 mov DL,AL ; 835 and AL,0xC3 ; 836 or AL,0x0C ; // round to 0 837 mov byte ptr cw+1,AL ; 838 fldcw cw ; 839 frndint ; 840 mov byte ptr cw+1,DL ; 841 fldcw cw ; 842 } 843 } 844 } 845 else 846 { 847 return core.stdc.math.truncl(x); 848 } 849 } 850 851 /// 852 @safe pure unittest 853 { 854 assert(trunc(0.01) == 0); 855 assert(trunc(0.49) == 0); 856 assert(trunc(0.5) == 0); 857 assert(trunc(1.5) == 1); 858 } 859 860 /***************************************** 861 * Returns x rounded to a long value using the current rounding mode. 862 * If the integer value of x is 863 * greater than long.max, the result is 864 * indeterminate. 865 */ 866 pragma(inline, true) 867 long rndtol(real x) @nogc @safe pure nothrow { return core.math.rndtol(x); } 868 //FIXME 869 ///ditto 870 pragma(inline, true) 871 long rndtol(double x) @safe pure nothrow @nogc { return rndtol(cast(real) x); } 872 //FIXME 873 ///ditto 874 pragma(inline, true) 875 long rndtol(float x) @safe pure nothrow @nogc { return rndtol(cast(real) x); } 876 877 /// 878 @safe unittest 879 { 880 assert(rndtol(1.0) == 1L); 881 assert(rndtol(1.2) == 1L); 882 assert(rndtol(1.7) == 2L); 883 assert(rndtol(1.0001) == 1L); 884 } 885 886 @safe unittest 887 { 888 long function(real) prndtol = &rndtol; 889 assert(prndtol != null); 890 } 891 892 // Helper for floor/ceil 893 T floorImpl(T)(const T x) @trusted pure nothrow @nogc 894 { 895 import std.math : floatTraits, RealFormat; 896 897 alias F = floatTraits!(T); 898 // Take care not to trigger library calls from the compiler, 899 // while ensuring that we don't get defeated by some optimizers. 900 union floatBits 901 { 902 T rv; 903 ushort[T.sizeof/2] vu; 904 905 // Other kinds of extractors for real formats. 906 static if (F.realFormat == RealFormat.ieeeSingle) 907 uint vi; 908 else static if (F.realFormat == RealFormat.ieeeDouble) 909 ulong vi; 910 } 911 floatBits y = void; 912 y.rv = x; 913 914 // Find the exponent (power of 2) 915 // Do this by shifting the raw value so that the exponent lies in the low bits, 916 // then mask out the sign bit, and subtract the bias. 917 static if (F.realFormat == RealFormat.ieeeSingle) 918 { 919 int exp = ((y.vi >> (T.mant_dig - 1)) & 0xff) - 0x7f; 920 enum mantissa_mask = F.MANTISSAMASK_INT; 921 enum sign_shift = 31; 922 } 923 else static if (F.realFormat == RealFormat.ieeeDouble) 924 { 925 long exp = ((y.vi >> (T.mant_dig - 1)) & 0x7ff) - 0x3ff; 926 enum mantissa_mask = F.MANTISSAMASK_LONG; 927 enum sign_shift = 63; 928 } 929 else static if (F.realFormat == RealFormat.ieeeExtended || 930 F.realFormat == RealFormat.ieeeExtended53) 931 { 932 int exp = (y.vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff; 933 934 version (LittleEndian) 935 int pos = 0; 936 else 937 int pos = 4; 938 } 939 else static if (F.realFormat == RealFormat.ieeeQuadruple) 940 { 941 int exp = (y.vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff; 942 943 version (LittleEndian) 944 int pos = 0; 945 else 946 int pos = 7; 947 } 948 else 949 static assert(false, "Not implemented for this architecture"); 950 951 if (exp < 0) 952 { 953 if (x < 0.0) 954 return -1.0; 955 else 956 return 0.0; 957 } 958 959 static if (F.realFormat == RealFormat.ieeeSingle || 960 F.realFormat == RealFormat.ieeeDouble) 961 { 962 if (exp < (T.mant_dig - 1)) 963 { 964 // Clear all bits representing the fraction part. 965 // Note: the fraction mask represents the floating point number 0.999999... 966 // i.e: `2.0 ^^ (exp - T.mant_dig + 1) * (fraction_mask + 1) == 1.0` 967 const fraction_mask = mantissa_mask >> exp; 968 969 if ((y.vi & fraction_mask) != 0) 970 { 971 // If 'x' is negative, then first substract (1.0 - T.epsilon) from the value. 972 if (y.vi >> sign_shift) 973 y.vi += fraction_mask; 974 y.vi &= ~fraction_mask; 975 } 976 } 977 } 978 else 979 { 980 static if (F.realFormat == RealFormat.ieeeExtended53) 981 exp = (T.mant_dig + 11 - 1) - exp; // mant_dig is really 64 982 else 983 exp = (T.mant_dig - 1) - exp; 984 985 // Zero 16 bits at a time. 986 while (exp >= 16) 987 { 988 version (LittleEndian) 989 y.vu[pos++] = 0; 990 else 991 y.vu[pos--] = 0; 992 exp -= 16; 993 } 994 995 // Clear the remaining bits. 996 if (exp > 0) 997 y.vu[pos] &= 0xffff ^ ((1 << exp) - 1); 998 999 if ((x < 0.0) && (x != y.rv)) 1000 y.rv -= 1.0; 1001 } 1002 1003 return y.rv; 1004 }