1 /*
2  * ====================================================
3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4  *
5  * Developed at SunPro, a Sun Microsystems, Inc. business.
6  * Permission to use, copy, modify, and distribute this
7  * software is freely granted, provided that this notice
8  * is preserved.
9  * ====================================================
10  */
11 /*
12 FUNCTION
13 <<round>>, <<roundf>>---round to integer, to nearest
14 INDEX
15 	round
16 INDEX
17 	roundf
18 
19 SYNOPSIS
20 	#include <math.h>
21 	double round(double <[x]>);
22 	float roundf(float <[x]>);
23 
24 DESCRIPTION
25 	The <<round>> functions round their argument to the nearest integer
26 	value in floating-point format, rounding halfway cases away from zero,
27 	regardless of the current rounding direction.  (While the "inexact"
28 	floating-point exception behavior is unspecified by the C standard, the
29 	<<round>> functions are written so that "inexact" is not raised if the
30 	result does not equal the argument, which behavior is as recommended by
31 	IEEE 754 for its related functions.)
32 
33 RETURNS
34 <[x]> rounded to an integral value.
35 
36 PORTABILITY
37 ANSI C, POSIX
38 
39 SEEALSO
40 <<nearbyint>>, <<rint>>
41 
42 */
43 
44 #include "fdlibm.h"
45 
46 #ifdef _NEED_FLOAT64
47 
48 __float64
round64(__float64 x)49 round64(__float64 x)
50 {
51   /* Most significant word, least significant word. */
52   __int32_t msw, exponent_less_1023;
53   __uint32_t lsw;
54 
55   EXTRACT_WORDS(msw, lsw, x);
56 
57   /* Extract exponent field. */
58   exponent_less_1023 = ((msw & 0x7ff00000) >> 20) - 1023;
59 
60   if (exponent_less_1023 < 20)
61     {
62       if (exponent_less_1023 < 0)
63         {
64           msw &= 0x80000000;
65           if (exponent_less_1023 == -1)
66             /* Result is +1.0 or -1.0. */
67             msw |= ((__int32_t)1023 << 20);
68           lsw = 0;
69         }
70       else
71         {
72           __uint32_t exponent_mask = 0x000fffff >> exponent_less_1023;
73           if ((msw & exponent_mask) == 0 && lsw == 0)
74             /* x in an integral value. */
75             return x;
76 
77           msw += 0x00080000 >> exponent_less_1023;
78           msw &= ~exponent_mask;
79           lsw = 0;
80         }
81     }
82   else if (exponent_less_1023 > 51)
83     {
84       if (exponent_less_1023 == 1024)
85         /* x is NaN or infinite. */
86         return x + x;
87       else
88         return x;
89     }
90   else
91     {
92       __uint32_t exponent_mask = 0xffffffff >> (exponent_less_1023 - 20);
93       __uint32_t tmp;
94 
95       if ((lsw & exponent_mask) == 0)
96         /* x is an integral value. */
97         return x;
98 
99       tmp = lsw + (1 << (51 - exponent_less_1023));
100       if (tmp < lsw)
101         msw += 1;
102       lsw = tmp;
103 
104       lsw &= ~exponent_mask;
105     }
106   INSERT_WORDS(x, msw, lsw);
107 
108   return x;
109 }
110 
111 _MATH_ALIAS_d_d(round)
112 
113 #endif /* _NEED_FLOAT64 */
114