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