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 | * rint(x) |
14 | * Return x rounded to integral value according to the prevailing |
15 | * rounding mode. |
16 | * Method: |
17 | * Using floating addition. |
18 | * Exception: |
19 | * Inexact flag raised if x not equal to rint(x). |
20 | */ |
21 | |
22 | #define NO_MATH_REDIRECT |
23 | #include <math.h> |
24 | #include <math_private.h> |
25 | #include <libm-alias-double.h> |
26 | #include <math-use-builtins.h> |
27 | |
28 | double |
29 | __rint (double x) |
30 | { |
31 | #if USE_RINT_BUILTIN |
32 | return __builtin_rint (x); |
33 | #else |
34 | /* Use generic implementation. */ |
35 | static const double |
36 | TWO52[2] = { |
37 | 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ |
38 | -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ |
39 | }; |
40 | int64_t i0, sx; |
41 | int32_t j0; |
42 | EXTRACT_WORDS64 (i0, x); |
43 | sx = (i0 >> 63) & 1; |
44 | j0 = ((i0 >> 52) & 0x7ff) - 0x3ff; |
45 | if (j0 < 52) |
46 | { |
47 | if (j0 < 0) |
48 | { |
49 | double w = TWO52[sx] + x; |
50 | double t = w - TWO52[sx]; |
51 | EXTRACT_WORDS64 (i0, t); |
52 | INSERT_WORDS64 (t, (i0 & UINT64_C (0x7fffffffffffffff)) |
53 | | (sx << 63)); |
54 | return t; |
55 | } |
56 | } |
57 | else |
58 | { |
59 | if (j0 == 0x400) |
60 | return x + x; /* inf or NaN */ |
61 | else |
62 | return x; /* x is integral */ |
63 | } |
64 | double w = TWO52[sx] + x; |
65 | return w - TWO52[sx]; |
66 | #endif /* ! USE_RINT_BUILTIN */ |
67 | } |
68 | #ifndef __rint |
69 | libm_alias_double (__rint, rint) |
70 | #endif |
71 |