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