1
2 /* @(#)s_rint.c 5.1 93/09/24 */
3 /*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13 /*
14 FUNCTION
15 <<rint>>, <<rintf>>---round to integer
16 INDEX
17 rint
18 INDEX
19 rintf
20
21 SYNOPSIS
22 #include <math.h>
23 double rint(double <[x]>);
24 float rintf(float <[x]>);
25
26 DESCRIPTION
27 The <<rint>> functions round their argument to an integer value in
28 floating-point format, using the current rounding direction. They
29 raise the "inexact" floating-point exception if the result differs
30 in value from the argument. See the <<nearbyint>> functions for the
31 same function with the "inexact" floating-point exception never being
32 raised. Newlib does not directly support floating-point exceptions.
33 The <<rint>> functions are written so that the "inexact" exception is
34 raised in hardware implementations that support it, even though Newlib
35 does not provide access.
36
37 RETURNS
38 <[x]> rounded to an integral value, using the current rounding direction.
39
40 PORTABILITY
41 ANSI C, POSIX
42
43 SEEALSO
44 <<nearbyint>>, <<round>>
45
46 */
47
48 /*
49 * rint(x)
50 * Return x rounded to integral value according to the prevailing
51 * rounding mode.
52 * Method:
53 * Using floating addition.
54 * Whenever a fraction is present, if the second or any following bit after
55 * the radix point is set, limit to the second radix point to avoid
56 * possible double rounding in the TWO52 +- steps (in case guard bits are
57 * used). Specifically, if have any, chop off bits past the 2nd place and
58 * set the second place.
59 * (e.g. 2.0625=0b10.0001 => 0b10.01=2.25;
60 * 2.3125=0b10.011 => 0b10.01=2.25;
61 * 1.5625= 0b1.1001 => 0b1.11=1.75;
62 * 1.9375= 0b1.1111 => 0b1.11=1.75.
63 * Pseudo-code: if(x.frac & ~0b0.10) x.frac = (x.frac & 0b0.11) | 0b0.01;).
64 * Exception:
65 * Inexact flag raised if x not equal to rint(x).
66 */
67
68 #include "fdlibm.h"
69
70 #ifdef _NEED_FLOAT64
71
72 #ifdef __STDC__
73 static const __float64
74 #else
75 static __float64
76 #endif
77 TWO52[2]={
78 _F_64(4.50359962737049600000e+15), /* 0x43300000, 0x00000000 */
79 _F_64(-4.50359962737049600000e+15), /* 0xC3300000, 0x00000000 */
80 };
81
82 __float64
rint64(__float64 x)83 rint64(__float64 x)
84 {
85 __int32_t i0,j0,sx;
86 __uint32_t i,i1;
87 __float64 t;
88 volatile __float64 w;
89 EXTRACT_WORDS(i0,i1,x);
90 sx = (i0>>31)&1; /* sign */
91 j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent */
92 if(j0<20) { /* no integral bits in LS part */
93 if(j0<0) { /* x is fractional or 0 */
94 if(((i0&0x7fffffff)|i1)==0) return x; /* x == 0 */
95 i1 |= (i0&0x0fffff);
96 i0 &= 0xfffe0000;
97 i0 |= ((i1|-i1)>>12)&0x80000;
98 SET_HIGH_WORD(x,i0);
99 w = TWO52[sx]+x;
100 t = w-TWO52[sx];
101 GET_HIGH_WORD(i0,t);
102 SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
103 return t;
104 } else { /* x has integer and maybe fraction */
105 i = (0x000fffff)>>j0;
106 if(((i0&i)|i1)==0) return x; /* x is integral */
107 i>>=1;
108 if(((i0&i)|i1)!=0) {
109 /* 2nd or any later bit after radix is set */
110 if(j0==19) i1 = 0x80000000; else i1 = 0;
111 i0 = (i0&(~i))|((0x40000)>>j0);
112 }
113 }
114 } else if (j0>51) {
115 /*
116 * Use barrier to avoid overflow on clang which would
117 * otherwise always do this add on arm and use a
118 * conditional move instead of a branch for the if
119 */
120 if (j0 == 0x400)
121 return opt_barrier_double(x+x);
122 return x;
123 } else {
124 i = ((__uint32_t)(0xffffffff))>>(j0-20);
125 if((i1&i)==0) return x; /* x is integral */
126 i>>=1;
127 if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
128 }
129 INSERT_WORDS(x,i0,i1);
130 w = TWO52[sx]+x;
131 return w-TWO52[sx];
132 }
133
134 _MATH_ALIAS_d_d(rint)
135
136 #endif /* _NEED_FLOAT64 */
137
