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 | #ifndef _DOUBLE_IS_32BITS |
---|
71 | |
---|
72 | #ifdef __STDC__ |
---|
73 | static const double |
---|
74 | #else |
---|
75 | static double |
---|
76 | #endif |
---|
77 | TWO52[2]={ |
---|
78 | 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ |
---|
79 | -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ |
---|
80 | }; |
---|
81 | |
---|
82 | #ifdef __STDC__ |
---|
83 | double rint(double x) |
---|
84 | #else |
---|
85 | double rint(x) |
---|
86 | double x; |
---|
87 | #endif |
---|
88 | { |
---|
89 | __int32_t i0,j0,sx; |
---|
90 | __uint32_t i,i1; |
---|
91 | double t; |
---|
92 | volatile double w; |
---|
93 | EXTRACT_WORDS(i0,i1,x); |
---|
94 | sx = (i0>>31)&1; /* sign */ |
---|
95 | j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent */ |
---|
96 | if(j0<20) { /* no integral bits in LS part */ |
---|
97 | if(j0<0) { /* x is fractional or 0 */ |
---|
98 | if(((i0&0x7fffffff)|i1)==0) return x; /* x == 0 */ |
---|
99 | i1 |= (i0&0x0fffff); |
---|
100 | i0 &= 0xfffe0000; |
---|
101 | i0 |= ((i1|-i1)>>12)&0x80000; |
---|
102 | SET_HIGH_WORD(x,i0); |
---|
103 | w = TWO52[sx]+x; |
---|
104 | t = w-TWO52[sx]; |
---|
105 | GET_HIGH_WORD(i0,t); |
---|
106 | SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31)); |
---|
107 | return t; |
---|
108 | } else { /* x has integer and maybe fraction */ |
---|
109 | i = (0x000fffff)>>j0; |
---|
110 | if(((i0&i)|i1)==0) return x; /* x is integral */ |
---|
111 | i>>=1; |
---|
112 | if(((i0&i)|i1)!=0) { |
---|
113 | /* 2nd or any later bit after radix is set */ |
---|
114 | if(j0==19) i1 = 0x80000000; else i1 = 0; |
---|
115 | i0 = (i0&(~i))|((0x40000)>>j0); |
---|
116 | } |
---|
117 | } |
---|
118 | } else if (j0>51) { |
---|
119 | if(j0==0x400) return x+x; /* inf or NaN */ |
---|
120 | else return x; /* x is integral */ |
---|
121 | } else { |
---|
122 | i = ((__uint32_t)(0xffffffff))>>(j0-20); |
---|
123 | if((i1&i)==0) return x; /* x is integral */ |
---|
124 | i>>=1; |
---|
125 | if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); |
---|
126 | } |
---|
127 | INSERT_WORDS(x,i0,i1); |
---|
128 | w = TWO52[sx]+x; |
---|
129 | return w-TWO52[sx]; |
---|
130 | } |
---|
131 | |
---|
132 | #endif /* _DOUBLE_IS_32BITS */ |
---|