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