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 | /* |
---|
13 | * Modified for ALMOS-MKH OS at UPMC, France, August 2018. (Alain Greiner) |
---|
14 | */ |
---|
15 | |
---|
16 | /* |
---|
17 | * rint(x) |
---|
18 | * Return x rounded to integral value according to the prevailing |
---|
19 | * rounding mode. |
---|
20 | * Method: |
---|
21 | * Using floating addition. |
---|
22 | * Exception: |
---|
23 | * Inexact flag raised if x not equal to rint(x). |
---|
24 | */ |
---|
25 | |
---|
26 | #include "math.h" |
---|
27 | #include "math_private.h" |
---|
28 | |
---|
29 | static const double |
---|
30 | TWO52[2]={ |
---|
31 | 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ |
---|
32 | -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ |
---|
33 | }; |
---|
34 | |
---|
35 | double rint(double x) |
---|
36 | { |
---|
37 | int32_t i0, j0, sx; |
---|
38 | uint32_t i,i1; |
---|
39 | double t; |
---|
40 | /* We use w = x + 2^52; t = w - 2^52; trick to round x to integer. |
---|
41 | * This trick requires that compiler does not optimize it |
---|
42 | * by keeping intermediate result w in a register wider than double. |
---|
43 | * Declaring w volatile assures that value gets truncated to double |
---|
44 | * (unfortunately, it also forces store+load): |
---|
45 | */ |
---|
46 | volatile double w; |
---|
47 | |
---|
48 | EXTRACT_WORDS(i0,i1,x); |
---|
49 | /* Unbiased exponent */ |
---|
50 | j0 = ((((uint32_t)i0) >> 20)&0x7ff)-0x3ff; |
---|
51 | |
---|
52 | if (j0 > 51) { |
---|
53 | //Why bother? Just returning x works too |
---|
54 | //if (j0 == 0x400) /* inf or NaN */ |
---|
55 | // return x+x; |
---|
56 | return x; /* x is integral */ |
---|
57 | } |
---|
58 | |
---|
59 | /* Sign */ |
---|
60 | sx = ((uint32_t)i0) >> 31; |
---|
61 | |
---|
62 | if (j0<20) { |
---|
63 | if (j0<0) { /* |x| < 1 */ |
---|
64 | if (((i0&0x7fffffff)|i1)==0) return x; |
---|
65 | i1 |= (i0&0x0fffff); |
---|
66 | i0 &= 0xfffe0000; |
---|
67 | i0 |= ((i1|-i1)>>12)&0x80000; |
---|
68 | SET_HIGH_WORD(x,i0); |
---|
69 | w = TWO52[sx]+x; |
---|
70 | t = w-TWO52[sx]; |
---|
71 | GET_HIGH_WORD(i0,t); |
---|
72 | SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31)); |
---|
73 | return t; |
---|
74 | } else { |
---|
75 | i = (0x000fffff)>>j0; |
---|
76 | if (((i0&i)|i1)==0) return x; /* x is integral */ |
---|
77 | i>>=1; |
---|
78 | if (((i0&i)|i1)!=0) { |
---|
79 | if (j0==19) i1 = 0x40000000; |
---|
80 | else i0 = (i0&(~i))|((0x20000)>>j0); |
---|
81 | } |
---|
82 | } |
---|
83 | } else { |
---|
84 | i = ((uint32_t)(0xffffffff))>>(j0-20); |
---|
85 | if ((i1&i)==0) return x; /* x is integral */ |
---|
86 | i>>=1; |
---|
87 | if ((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); |
---|
88 | } |
---|
89 | INSERT_WORDS(x,i0,i1); |
---|
90 | w = TWO52[sx]+x; |
---|
91 | return w-TWO52[sx]; |
---|
92 | } |
---|