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 /* NOTE: This file has been patched from the original DMD distribution to
22 * work with the GDC compiler.
23 */
24 module std.math.rounding;
25
26 static import core.math;
27 static import core.stdc.math;
28
29 import std.traits : isFloatingPoint, isIntegral, Unqual;
30
31 version (D_InlineAsm_X86) version = InlineAsm_X86_Any;
32 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any;
33
34 version (InlineAsm_X86_Any) version = InlineAsm_X87;
version(InlineAsm_X87)35 version (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 */
ceil(real x)45 real 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
ceil(double x)121 double 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
ceil(float x)153 float 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 */
floor(real x)188 real 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
floor(double x)262 double 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
floor(float x)292 float 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 */
334 Unqual!F quantize(alias rfunc = rint, F)(const F val, const F unit)
335 if (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 */
380 Unqual!F quantize(real base, alias rfunc = rint, F, E)(const F val, const E exp)
381 if (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
390 Unqual!F quantize(real base, long exp = 1, alias rfunc = rint, F)(const F val)
391 if (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 */
442 pragma(inline, true)
443 real 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 */
473 pragma(inline, true)
474 real rint(real x) @safe pure nothrow @nogc
475 {
476 return core.math.rint(x);
477 }
478 ///ditto
479 pragma(inline, true)
480 double rint(double x) @safe pure nothrow @nogc
481 {
482 return core.math.rint(x);
483 }
484 ///ditto
485 pragma(inline, true)
486 float 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 */
526 long 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
717 static 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 */
736 auto 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
766 version (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 */
782 long 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 */
807 real 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 */
870 pragma(inline, true)
871 long rndtol(real x) @nogc @safe pure nothrow { return core.math.rndtol(x); }
872 //FIXME
873 ///ditto
874 pragma(inline, true)
875 long rndtol(double x) @safe pure nothrow @nogc { return rndtol(cast(real) x); }
876 //FIXME
877 ///ditto
878 pragma(inline, true)
879 long 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
897 T 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