]>
Commit | Line | Data |
---|---|---|
f964490f RM |
1 | #ifndef _MATH_PRIVATE_H_ |
2 | #error "Never use <math_ldbl.h> directly; include <math_private.h> instead." | |
3 | #endif | |
4 | ||
5 | #include <sysdeps/ieee754/ldbl-128/math_ldbl.h> | |
5c68d401 | 6 | #include <ieee754.h> |
e054f494 | 7 | #include <stdint.h> |
9c84384c | 8 | |
5c68d401 | 9 | static inline void |
34ae0b32 | 10 | ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x) |
5c68d401 RM |
11 | { |
12 | /* We have 105 bits of mantissa plus one implicit digit. Since | |
13 | 106 bits are representable we use the first implicit digit for | |
14 | the number before the decimal point and the second implicit bit | |
15 | as bit 53 of the mantissa. */ | |
34ae0b32 | 16 | uint64_t hi, lo; |
5c68d401 | 17 | int ediff; |
9605ca6c AM |
18 | union ibm_extended_long_double u; |
19 | u.ld = x; | |
20 | *exp = u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS; | |
f964490f | 21 | |
9605ca6c AM |
22 | lo = ((uint64_t) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1; |
23 | hi = ((uint64_t) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1; | |
5c68d401 RM |
24 | /* If the lower double is not a denomal or zero then set the hidden |
25 | 53rd bit. */ | |
9605ca6c | 26 | if (u.d[1].ieee.exponent > 0x001) |
5c68d401 RM |
27 | { |
28 | lo |= (1ULL << 52); | |
29 | lo = lo << 7; /* pre-shift lo to match ieee854. */ | |
30 | /* The lower double is normalized separately from the upper. We | |
31 | may need to adjust the lower manitissa to reflect this. */ | |
9605ca6c | 32 | ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent; |
5c68d401 RM |
33 | if (ediff > 53) |
34 | lo = lo >> (ediff-53); | |
34ae0b32 | 35 | hi |= (1ULL << 52); |
5c68d401 | 36 | } |
9c84384c | 37 | |
9605ca6c AM |
38 | if ((u.d[0].ieee.negative != u.d[1].ieee.negative) |
39 | && ((u.d[1].ieee.exponent != 0) && (lo != 0LL))) | |
5c68d401 RM |
40 | { |
41 | hi--; | |
42 | lo = (1ULL << 60) - lo; | |
43 | if (hi < (1ULL << 52)) | |
44 | { | |
45 | /* we have a borrow from the hidden bit, so shift left 1. */ | |
46 | hi = (hi << 1) | (lo >> 59); | |
47 | lo = 0xfffffffffffffffLL & (lo << 1); | |
48 | *exp = *exp - 1; | |
49 | } | |
50 | } | |
51 | *lo64 = (hi << 60) | lo; | |
52 | *hi64 = hi >> 4; | |
53 | } | |
f964490f | 54 | |
5c68d401 RM |
55 | static inline long double |
56 | ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64) | |
57 | { | |
58 | union ibm_extended_long_double u; | |
59 | unsigned long hidden2, lzcount; | |
60 | unsigned long long hi, lo; | |
61 | ||
9605ca6c AM |
62 | u.d[0].ieee.negative = sign; |
63 | u.d[1].ieee.negative = sign; | |
64 | u.d[0].ieee.exponent = exp + IEEE754_DOUBLE_BIAS; | |
65 | u.d[1].ieee.exponent = exp-53 + IEEE754_DOUBLE_BIAS; | |
5c68d401 | 66 | /* Expect 113 bits (112 bits + hidden) right justified in two longs. |
9c84384c | 67 | The low order 53 bits (52 + hidden) go into the lower double */ |
5c68d401 RM |
68 | lo = (lo64 >> 7)& ((1ULL << 53) - 1); |
69 | hidden2 = (lo64 >> 59) & 1ULL; | |
70 | /* The high order 53 bits (52 + hidden) go into the upper double */ | |
71 | hi = (lo64 >> 60) & ((1ULL << 11) - 1); | |
72 | hi |= (hi64 << 4); | |
73 | ||
74 | if (lo != 0LL) | |
75 | { | |
76 | /* hidden2 bit of low double controls rounding of the high double. | |
77 | If hidden2 is '1' then round up hi and adjust lo (2nd mantissa) | |
78 | plus change the sign of the low double to compensate. */ | |
79 | if (hidden2) | |
80 | { | |
81 | hi++; | |
9605ca6c | 82 | u.d[1].ieee.negative = !sign; |
5c68d401 RM |
83 | lo = (1ULL << 53) - lo; |
84 | } | |
85 | /* The hidden bit of the lo mantissa is zero so we need to | |
86 | normalize the it for the low double. Shift it left until the | |
9c84384c | 87 | hidden bit is '1' then adjust the 2nd exponent accordingly. */ |
5c68d401 RM |
88 | |
89 | if (sizeof (lo) == sizeof (long)) | |
90 | lzcount = __builtin_clzl (lo); | |
91 | else if ((lo >> 32) != 0) | |
92 | lzcount = __builtin_clzl ((long) (lo >> 32)); | |
93 | else | |
94 | lzcount = __builtin_clzl ((long) lo) + 32; | |
95 | lzcount = lzcount - 11; | |
96 | if (lzcount > 0) | |
97 | { | |
9605ca6c | 98 | int expnt2 = u.d[1].ieee.exponent - lzcount; |
5c68d401 RM |
99 | if (expnt2 >= 1) |
100 | { | |
101 | /* Not denormal. Normalize and set low exponent. */ | |
102 | lo = lo << lzcount; | |
9605ca6c | 103 | u.d[1].ieee.exponent = expnt2; |
5c68d401 RM |
104 | } |
105 | else | |
106 | { | |
107 | /* Is denormal. */ | |
108 | lo = lo << (lzcount + expnt2); | |
9605ca6c | 109 | u.d[1].ieee.exponent = 0; |
5c68d401 RM |
110 | } |
111 | } | |
112 | } | |
113 | else | |
114 | { | |
9605ca6c AM |
115 | u.d[1].ieee.negative = 0; |
116 | u.d[1].ieee.exponent = 0; | |
5c68d401 RM |
117 | } |
118 | ||
9605ca6c AM |
119 | u.d[1].ieee.mantissa1 = lo & ((1ULL << 32) - 1); |
120 | u.d[1].ieee.mantissa0 = (lo >> 32) & ((1ULL << 20) - 1); | |
121 | u.d[0].ieee.mantissa1 = hi & ((1ULL << 32) - 1); | |
122 | u.d[0].ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1); | |
123 | return u.ld; | |
5c68d401 | 124 | } |
9c84384c | 125 | |
5c68d401 RM |
126 | /* Handy utility functions to pack/unpack/cononicalize and find the nearbyint |
127 | of long double implemented as double double. */ | |
128 | static inline long double | |
edf66e57 | 129 | default_ldbl_pack (double a, double aa) |
5c68d401 RM |
130 | { |
131 | union ibm_extended_long_double u; | |
9605ca6c AM |
132 | u.d[0].d = a; |
133 | u.d[1].d = aa; | |
134 | return u.ld; | |
5c68d401 RM |
135 | } |
136 | ||
137 | static inline void | |
edf66e57 | 138 | default_ldbl_unpack (long double l, double *a, double *aa) |
5c68d401 RM |
139 | { |
140 | union ibm_extended_long_double u; | |
9605ca6c AM |
141 | u.ld = l; |
142 | *a = u.d[0].d; | |
143 | *aa = u.d[1].d; | |
5c68d401 RM |
144 | } |
145 | ||
edf66e57 AZ |
146 | #ifndef ldbl_pack |
147 | # define ldbl_pack default_ldbl_pack | |
148 | #endif | |
149 | #ifndef ldbl_unpack | |
150 | # define ldbl_unpack default_ldbl_unpack | |
151 | #endif | |
5c68d401 RM |
152 | |
153 | /* Convert a finite long double to canonical form. | |
154 | Does not handle +/-Inf properly. */ | |
155 | static inline void | |
156 | ldbl_canonicalize (double *a, double *aa) | |
157 | { | |
158 | double xh, xl; | |
159 | ||
160 | xh = *a + *aa; | |
161 | xl = (*a - xh) + *aa; | |
162 | *a = xh; | |
163 | *aa = xl; | |
164 | } | |
165 | ||
166 | /* Simple inline nearbyint (double) function . | |
167 | Only works in the default rounding mode | |
168 | but is useful in long double rounding functions. */ | |
169 | static inline double | |
170 | ldbl_nearbyint (double a) | |
171 | { | |
172 | double two52 = 0x10000000000000LL; | |
173 | ||
174 | if (__builtin_expect ((__builtin_fabs (a) < two52), 1)) | |
175 | { | |
176 | if (__builtin_expect ((a > 0.0), 1)) | |
177 | { | |
178 | a += two52; | |
179 | a -= two52; | |
180 | } | |
181 | else if (__builtin_expect ((a < 0.0), 1)) | |
182 | { | |
183 | a = two52 - a; | |
184 | a = -(a - two52); | |
185 | } | |
186 | } | |
187 | return a; | |
188 | } |