|
0
|
1 /* Copyright (C) 2007, 2009 Free Software Foundation, Inc.
|
|
|
2
|
|
|
3 This file is part of GCC.
|
|
|
4
|
|
|
5 GCC is free software; you can redistribute it and/or modify it under
|
|
|
6 the terms of the GNU General Public License as published by the Free
|
|
|
7 Software Foundation; either version 3, or (at your option) any later
|
|
|
8 version.
|
|
|
9
|
|
|
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
|
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
|
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
|
13 for more details.
|
|
|
14
|
|
|
15 Under Section 7 of GPL version 3, you are granted additional
|
|
|
16 permissions described in the GCC Runtime Library Exception, version
|
|
|
17 3.1, as published by the Free Software Foundation.
|
|
|
18
|
|
|
19 You should have received a copy of the GNU General Public License and
|
|
|
20 a copy of the GCC Runtime Library Exception along with this program;
|
|
|
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
|
|
22 <http://www.gnu.org/licenses/>. */
|
|
|
23
|
|
|
24 #define BID_128RES
|
|
|
25 #include "bid_internal.h"
|
|
|
26
|
|
|
27 /*****************************************************************************
|
|
|
28 * BID128 nextup
|
|
|
29 ****************************************************************************/
|
|
|
30
|
|
|
31 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
32 void
|
|
|
33 bid128_nextup (UINT128 * pres,
|
|
|
34 UINT128 *
|
|
|
35 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
|
|
36 UINT128 x = *px;
|
|
|
37 #else
|
|
|
38 UINT128
|
|
|
39 bid128_nextup (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
|
|
40 _EXC_INFO_PARAM) {
|
|
|
41 #endif
|
|
|
42
|
|
|
43 UINT128 res;
|
|
|
44 UINT64 x_sign;
|
|
|
45 UINT64 x_exp;
|
|
|
46 int exp;
|
|
|
47 BID_UI64DOUBLE tmp1;
|
|
|
48 int x_nr_bits;
|
|
|
49 int q1, ind;
|
|
|
50 UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
|
|
|
51
|
|
|
52 BID_SWAP128 (x);
|
|
|
53 // unpack the argument
|
|
|
54 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
|
|
|
55 C1.w[1] = x.w[1] & MASK_COEFF;
|
|
|
56 C1.w[0] = x.w[0];
|
|
|
57
|
|
|
58 // check for NaN or Infinity
|
|
|
59 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
|
|
|
60 // x is special
|
|
|
61 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
|
|
|
62 // if x = NaN, then res = Q (x)
|
|
|
63 // check first for non-canonical NaN payload
|
|
|
64 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
|
|
|
65 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
|
|
|
66 && (x.w[0] > 0x38c15b09ffffffffull))) {
|
|
|
67 x.w[1] = x.w[1] & 0xffffc00000000000ull;
|
|
|
68 x.w[0] = 0x0ull;
|
|
|
69 }
|
|
|
70 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
|
|
|
71 // set invalid flag
|
|
|
72 *pfpsf |= INVALID_EXCEPTION;
|
|
|
73 // return quiet (x)
|
|
|
74 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
|
|
|
75 res.w[0] = x.w[0];
|
|
|
76 } else { // x is QNaN
|
|
|
77 // return x
|
|
|
78 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
|
|
|
79 res.w[0] = x.w[0];
|
|
|
80 }
|
|
|
81 } else { // x is not NaN, so it must be infinity
|
|
|
82 if (!x_sign) { // x is +inf
|
|
|
83 res.w[1] = 0x7800000000000000ull; // +inf
|
|
|
84 res.w[0] = 0x0000000000000000ull;
|
|
|
85 } else { // x is -inf
|
|
|
86 res.w[1] = 0xdfffed09bead87c0ull; // -MAXFP = -999...99 * 10^emax
|
|
|
87 res.w[0] = 0x378d8e63ffffffffull;
|
|
|
88 }
|
|
|
89 }
|
|
|
90 BID_RETURN (res);
|
|
|
91 }
|
|
|
92 // check for non-canonical values (treated as zero)
|
|
|
93 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
|
|
|
94 // non-canonical
|
|
|
95 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
|
|
|
96 C1.w[1] = 0; // significand high
|
|
|
97 C1.w[0] = 0; // significand low
|
|
|
98 } else { // G0_G1 != 11
|
|
|
99 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
|
|
|
100 if (C1.w[1] > 0x0001ed09bead87c0ull ||
|
|
|
101 (C1.w[1] == 0x0001ed09bead87c0ull
|
|
|
102 && C1.w[0] > 0x378d8e63ffffffffull)) {
|
|
|
103 // x is non-canonical if coefficient is larger than 10^34 -1
|
|
|
104 C1.w[1] = 0;
|
|
|
105 C1.w[0] = 0;
|
|
|
106 } else { // canonical
|
|
|
107 ;
|
|
|
108 }
|
|
|
109 }
|
|
|
110
|
|
|
111 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
|
|
|
112 // x is +/-0
|
|
|
113 res.w[1] = 0x0000000000000000ull; // +1 * 10^emin
|
|
|
114 res.w[0] = 0x0000000000000001ull;
|
|
|
115 } else { // x is not special and is not zero
|
|
|
116 if (x.w[1] == 0x5fffed09bead87c0ull
|
|
|
117 && x.w[0] == 0x378d8e63ffffffffull) {
|
|
|
118 // x = +MAXFP = 999...99 * 10^emax
|
|
|
119 res.w[1] = 0x7800000000000000ull; // +inf
|
|
|
120 res.w[0] = 0x0000000000000000ull;
|
|
|
121 } else if (x.w[1] == 0x8000000000000000ull
|
|
|
122 && x.w[0] == 0x0000000000000001ull) {
|
|
|
123 // x = -MINFP = 1...99 * 10^emin
|
|
|
124 res.w[1] = 0x8000000000000000ull; // -0
|
|
|
125 res.w[0] = 0x0000000000000000ull;
|
|
|
126 } else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
|
|
127 // can add/subtract 1 ulp to the significand
|
|
|
128
|
|
|
129 // Note: we could check here if x >= 10^34 to speed up the case q1 = 34
|
|
|
130 // q1 = nr. of decimal digits in x
|
|
|
131 // determine first the nr. of bits in x
|
|
|
132 if (C1.w[1] == 0) {
|
|
|
133 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
|
|
|
134 // split the 64-bit value in two 32-bit halves to avoid rnd errors
|
|
|
135 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
|
|
|
136 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
|
|
|
137 x_nr_bits =
|
|
|
138 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
|
|
|
139 0x3ff);
|
|
|
140 } else { // x < 2^32
|
|
|
141 tmp1.d = (double) (C1.w[0]); // exact conversion
|
|
|
142 x_nr_bits =
|
|
|
143 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
|
|
|
144 0x3ff);
|
|
|
145 }
|
|
|
146 } else { // if x < 2^53
|
|
|
147 tmp1.d = (double) C1.w[0]; // exact conversion
|
|
|
148 x_nr_bits =
|
|
|
149 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
|
150 }
|
|
|
151 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
|
|
|
152 tmp1.d = (double) C1.w[1]; // exact conversion
|
|
|
153 x_nr_bits =
|
|
|
154 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
|
155 }
|
|
|
156 q1 = nr_digits[x_nr_bits - 1].digits;
|
|
|
157 if (q1 == 0) {
|
|
|
158 q1 = nr_digits[x_nr_bits - 1].digits1;
|
|
|
159 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
|
|
|
160 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
|
|
|
161 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
|
|
|
162 q1++;
|
|
|
163 }
|
|
|
164 // if q1 < P34 then pad the significand with zeros
|
|
|
165 if (q1 < P34) {
|
|
|
166 exp = (x_exp >> 49) - 6176;
|
|
|
167 if (exp + 6176 > P34 - q1) {
|
|
|
168 ind = P34 - q1; // 1 <= ind <= P34 - 1
|
|
|
169 // pad with P34 - q1 zeros, until exponent = emin
|
|
|
170 // C1 = C1 * 10^ind
|
|
|
171 if (q1 <= 19) { // 64-bit C1
|
|
|
172 if (ind <= 19) { // 64-bit 10^ind and 64-bit C1
|
|
|
173 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
|
|
174 } else { // 128-bit 10^ind and 64-bit C1
|
|
|
175 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
|
|
176 }
|
|
|
177 } else { // C1 is (most likely) 128-bit
|
|
|
178 if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely)
|
|
|
179 __mul_128x64_to_128 (C1, ten2k64[ind], C1);
|
|
|
180 } else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
|
|
|
181 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
|
|
182 } else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
|
|
|
183 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
|
|
184 }
|
|
|
185 }
|
|
|
186 x_exp = x_exp - ((UINT64) ind << 49);
|
|
|
187 } else { // pad with zeros until the exponent reaches emin
|
|
|
188 ind = exp + 6176;
|
|
|
189 // C1 = C1 * 10^ind
|
|
|
190 if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
|
|
|
191 if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind
|
|
|
192 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
|
|
193 } else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
|
|
|
194 __mul_128x64_to_128 (C1, ten2k64[ind], C1);
|
|
|
195 }
|
|
|
196 } else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
|
|
|
197 // 64-bit C1, 128-bit 10^ind
|
|
|
198 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
|
|
199 }
|
|
|
200 x_exp = EXP_MIN;
|
|
|
201 }
|
|
|
202 }
|
|
|
203 if (!x_sign) { // x > 0
|
|
|
204 // add 1 ulp (add 1 to the significand)
|
|
|
205 C1.w[0]++;
|
|
|
206 if (C1.w[0] == 0)
|
|
|
207 C1.w[1]++;
|
|
|
208 if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34
|
|
|
209 C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33
|
|
|
210 C1.w[0] = 0x38c15b0a00000000ull;
|
|
|
211 x_exp = x_exp + EXP_P1;
|
|
|
212 }
|
|
|
213 } else { // x < 0
|
|
|
214 // subtract 1 ulp (subtract 1 from the significand)
|
|
|
215 C1.w[0]--;
|
|
|
216 if (C1.w[0] == 0xffffffffffffffffull)
|
|
|
217 C1.w[1]--;
|
|
|
218 if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1
|
|
|
219 C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1
|
|
|
220 C1.w[0] = 0x378d8e63ffffffffull;
|
|
|
221 x_exp = x_exp - EXP_P1;
|
|
|
222 }
|
|
|
223 }
|
|
|
224 // assemble the result
|
|
|
225 res.w[1] = x_sign | x_exp | C1.w[1];
|
|
|
226 res.w[0] = C1.w[0];
|
|
|
227 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
|
|
228 } // end x is not special and is not zero
|
|
|
229 BID_RETURN (res);
|
|
|
230 }
|
|
|
231
|
|
|
232 /*****************************************************************************
|
|
|
233 * BID128 nextdown
|
|
|
234 ****************************************************************************/
|
|
|
235
|
|
|
236 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
237 void
|
|
|
238 bid128_nextdown (UINT128 * pres,
|
|
|
239 UINT128 *
|
|
|
240 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
|
|
241 UINT128 x = *px;
|
|
|
242 #else
|
|
|
243 UINT128
|
|
|
244 bid128_nextdown (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
|
|
245 _EXC_INFO_PARAM) {
|
|
|
246 #endif
|
|
|
247
|
|
|
248 UINT128 res;
|
|
|
249 UINT64 x_sign;
|
|
|
250 UINT64 x_exp;
|
|
|
251 int exp;
|
|
|
252 BID_UI64DOUBLE tmp1;
|
|
|
253 int x_nr_bits;
|
|
|
254 int q1, ind;
|
|
|
255 UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
|
|
|
256
|
|
|
257 BID_SWAP128 (x);
|
|
|
258 // unpack the argument
|
|
|
259 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
|
|
|
260 C1.w[1] = x.w[1] & MASK_COEFF;
|
|
|
261 C1.w[0] = x.w[0];
|
|
|
262
|
|
|
263 // check for NaN or Infinity
|
|
|
264 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
|
|
|
265 // x is special
|
|
|
266 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
|
|
|
267 // if x = NaN, then res = Q (x)
|
|
|
268 // check first for non-canonical NaN payload
|
|
|
269 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
|
|
|
270 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
|
|
|
271 && (x.w[0] > 0x38c15b09ffffffffull))) {
|
|
|
272 x.w[1] = x.w[1] & 0xffffc00000000000ull;
|
|
|
273 x.w[0] = 0x0ull;
|
|
|
274 }
|
|
|
275 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
|
|
|
276 // set invalid flag
|
|
|
277 *pfpsf |= INVALID_EXCEPTION;
|
|
|
278 // return quiet (x)
|
|
|
279 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
|
|
|
280 res.w[0] = x.w[0];
|
|
|
281 } else { // x is QNaN
|
|
|
282 // return x
|
|
|
283 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
|
|
|
284 res.w[0] = x.w[0];
|
|
|
285 }
|
|
|
286 } else { // x is not NaN, so it must be infinity
|
|
|
287 if (!x_sign) { // x is +inf
|
|
|
288 res.w[1] = 0x5fffed09bead87c0ull; // +MAXFP = +999...99 * 10^emax
|
|
|
289 res.w[0] = 0x378d8e63ffffffffull;
|
|
|
290 } else { // x is -inf
|
|
|
291 res.w[1] = 0xf800000000000000ull; // -inf
|
|
|
292 res.w[0] = 0x0000000000000000ull;
|
|
|
293 }
|
|
|
294 }
|
|
|
295 BID_RETURN (res);
|
|
|
296 }
|
|
|
297 // check for non-canonical values (treated as zero)
|
|
|
298 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
|
|
|
299 // non-canonical
|
|
|
300 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
|
|
|
301 C1.w[1] = 0; // significand high
|
|
|
302 C1.w[0] = 0; // significand low
|
|
|
303 } else { // G0_G1 != 11
|
|
|
304 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
|
|
|
305 if (C1.w[1] > 0x0001ed09bead87c0ull ||
|
|
|
306 (C1.w[1] == 0x0001ed09bead87c0ull
|
|
|
307 && C1.w[0] > 0x378d8e63ffffffffull)) {
|
|
|
308 // x is non-canonical if coefficient is larger than 10^34 -1
|
|
|
309 C1.w[1] = 0;
|
|
|
310 C1.w[0] = 0;
|
|
|
311 } else { // canonical
|
|
|
312 ;
|
|
|
313 }
|
|
|
314 }
|
|
|
315
|
|
|
316 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
|
|
|
317 // x is +/-0
|
|
|
318 res.w[1] = 0x8000000000000000ull; // -1 * 10^emin
|
|
|
319 res.w[0] = 0x0000000000000001ull;
|
|
|
320 } else { // x is not special and is not zero
|
|
|
321 if (x.w[1] == 0xdfffed09bead87c0ull
|
|
|
322 && x.w[0] == 0x378d8e63ffffffffull) {
|
|
|
323 // x = -MAXFP = -999...99 * 10^emax
|
|
|
324 res.w[1] = 0xf800000000000000ull; // -inf
|
|
|
325 res.w[0] = 0x0000000000000000ull;
|
|
|
326 } else if (x.w[1] == 0x0ull && x.w[0] == 0x0000000000000001ull) { // +MINFP
|
|
|
327 res.w[1] = 0x0000000000000000ull; // +0
|
|
|
328 res.w[0] = 0x0000000000000000ull;
|
|
|
329 } else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
|
|
330 // can add/subtract 1 ulp to the significand
|
|
|
331
|
|
|
332 // Note: we could check here if x >= 10^34 to speed up the case q1 = 34
|
|
|
333 // q1 = nr. of decimal digits in x
|
|
|
334 // determine first the nr. of bits in x
|
|
|
335 if (C1.w[1] == 0) {
|
|
|
336 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
|
|
|
337 // split the 64-bit value in two 32-bit halves to avoid rnd errors
|
|
|
338 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
|
|
|
339 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
|
|
|
340 x_nr_bits =
|
|
|
341 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
|
|
|
342 0x3ff);
|
|
|
343 } else { // x < 2^32
|
|
|
344 tmp1.d = (double) (C1.w[0]); // exact conversion
|
|
|
345 x_nr_bits =
|
|
|
346 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
|
|
|
347 0x3ff);
|
|
|
348 }
|
|
|
349 } else { // if x < 2^53
|
|
|
350 tmp1.d = (double) C1.w[0]; // exact conversion
|
|
|
351 x_nr_bits =
|
|
|
352 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
|
353 }
|
|
|
354 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
|
|
|
355 tmp1.d = (double) C1.w[1]; // exact conversion
|
|
|
356 x_nr_bits =
|
|
|
357 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
|
358 }
|
|
|
359 q1 = nr_digits[x_nr_bits - 1].digits;
|
|
|
360 if (q1 == 0) {
|
|
|
361 q1 = nr_digits[x_nr_bits - 1].digits1;
|
|
|
362 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
|
|
|
363 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
|
|
|
364 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
|
|
|
365 q1++;
|
|
|
366 }
|
|
|
367 // if q1 < P then pad the significand with zeros
|
|
|
368 if (q1 < P34) {
|
|
|
369 exp = (x_exp >> 49) - 6176;
|
|
|
370 if (exp + 6176 > P34 - q1) {
|
|
|
371 ind = P34 - q1; // 1 <= ind <= P34 - 1
|
|
|
372 // pad with P34 - q1 zeros, until exponent = emin
|
|
|
373 // C1 = C1 * 10^ind
|
|
|
374 if (q1 <= 19) { // 64-bit C1
|
|
|
375 if (ind <= 19) { // 64-bit 10^ind and 64-bit C1
|
|
|
376 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
|
|
377 } else { // 128-bit 10^ind and 64-bit C1
|
|
|
378 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
|
|
379 }
|
|
|
380 } else { // C1 is (most likely) 128-bit
|
|
|
381 if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely)
|
|
|
382 __mul_128x64_to_128 (C1, ten2k64[ind], C1);
|
|
|
383 } else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
|
|
|
384 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
|
|
385 } else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
|
|
|
386 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
|
|
387 }
|
|
|
388 }
|
|
|
389 x_exp = x_exp - ((UINT64) ind << 49);
|
|
|
390 } else { // pad with zeros until the exponent reaches emin
|
|
|
391 ind = exp + 6176;
|
|
|
392 // C1 = C1 * 10^ind
|
|
|
393 if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
|
|
|
394 if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind
|
|
|
395 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
|
|
396 } else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
|
|
|
397 __mul_128x64_to_128 (C1, ten2k64[ind], C1);
|
|
|
398 }
|
|
|
399 } else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
|
|
|
400 // 64-bit C1, 128-bit 10^ind
|
|
|
401 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
|
|
402 }
|
|
|
403 x_exp = EXP_MIN;
|
|
|
404 }
|
|
|
405 }
|
|
|
406 if (x_sign) { // x < 0
|
|
|
407 // add 1 ulp (add 1 to the significand)
|
|
|
408 C1.w[0]++;
|
|
|
409 if (C1.w[0] == 0)
|
|
|
410 C1.w[1]++;
|
|
|
411 if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34
|
|
|
412 C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33
|
|
|
413 C1.w[0] = 0x38c15b0a00000000ull;
|
|
|
414 x_exp = x_exp + EXP_P1;
|
|
|
415 }
|
|
|
416 } else { // x > 0
|
|
|
417 // subtract 1 ulp (subtract 1 from the significand)
|
|
|
418 C1.w[0]--;
|
|
|
419 if (C1.w[0] == 0xffffffffffffffffull)
|
|
|
420 C1.w[1]--;
|
|
|
421 if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1
|
|
|
422 C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1
|
|
|
423 C1.w[0] = 0x378d8e63ffffffffull;
|
|
|
424 x_exp = x_exp - EXP_P1;
|
|
|
425 }
|
|
|
426 }
|
|
|
427 // assemble the result
|
|
|
428 res.w[1] = x_sign | x_exp | C1.w[1];
|
|
|
429 res.w[0] = C1.w[0];
|
|
|
430 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
|
|
431 } // end x is not special and is not zero
|
|
|
432 BID_RETURN (res);
|
|
|
433 }
|
|
|
434
|
|
|
435 /*****************************************************************************
|
|
|
436 * BID128 nextafter
|
|
|
437 ****************************************************************************/
|
|
|
438
|
|
|
439 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
440 void
|
|
|
441 bid128_nextafter (UINT128 * pres, UINT128 * px,
|
|
|
442 UINT128 *
|
|
|
443 py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)
|
|
|
444 {
|
|
|
445 UINT128 x = *px;
|
|
|
446 UINT128 y = *py;
|
|
|
447 UINT128 xnswp = *px;
|
|
|
448 UINT128 ynswp = *py;
|
|
|
449 #else
|
|
|
450 UINT128
|
|
|
451 bid128_nextafter (UINT128 x,
|
|
|
452 UINT128 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
|
|
453 _EXC_INFO_PARAM) {
|
|
|
454 UINT128 xnswp = x;
|
|
|
455 UINT128 ynswp = y;
|
|
|
456 #endif
|
|
|
457
|
|
|
458 UINT128 res;
|
|
|
459 UINT128 tmp1, tmp2, tmp3;
|
|
|
460 FPSC tmp_fpsf = 0; // dummy fpsf for calls to comparison functions
|
|
|
461 int res1, res2;
|
|
|
462 UINT64 x_exp;
|
|
|
463
|
|
|
464
|
|
|
465 BID_SWAP128 (x);
|
|
|
466 BID_SWAP128 (y);
|
|
|
467 // check for NaNs
|
|
|
468 if (((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL)
|
|
|
469 || ((y.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
|
|
|
470 // x is special or y is special
|
|
|
471 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
|
|
|
472 // if x = NaN, then res = Q (x)
|
|
|
473 // check first for non-canonical NaN payload
|
|
|
474 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
|
|
|
475 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
|
|
|
476 && (x.w[0] > 0x38c15b09ffffffffull))) {
|
|
|
477 x.w[1] = x.w[1] & 0xffffc00000000000ull;
|
|
|
478 x.w[0] = 0x0ull;
|
|
|
479 }
|
|
|
480 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
|
|
|
481 // set invalid flag
|
|
|
482 *pfpsf |= INVALID_EXCEPTION;
|
|
|
483 // return quiet (x)
|
|
|
484 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
|
|
|
485 res.w[0] = x.w[0];
|
|
|
486 } else { // x is QNaN
|
|
|
487 // return x
|
|
|
488 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
|
|
|
489 res.w[0] = x.w[0];
|
|
|
490 if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
|
|
|
491 // set invalid flag
|
|
|
492 *pfpsf |= INVALID_EXCEPTION;
|
|
|
493 }
|
|
|
494 }
|
|
|
495 BID_RETURN (res)
|
|
|
496 } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
|
|
|
497 // if x = NaN, then res = Q (x)
|
|
|
498 // check first for non-canonical NaN payload
|
|
|
499 if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
|
|
|
500 (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
|
|
|
501 && (y.w[0] > 0x38c15b09ffffffffull))) {
|
|
|
502 y.w[1] = y.w[1] & 0xffffc00000000000ull;
|
|
|
503 y.w[0] = 0x0ull;
|
|
|
504 }
|
|
|
505 if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
|
|
|
506 // set invalid flag
|
|
|
507 *pfpsf |= INVALID_EXCEPTION;
|
|
|
508 // return quiet (x)
|
|
|
509 res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
|
|
|
510 res.w[0] = y.w[0];
|
|
|
511 } else { // x is QNaN
|
|
|
512 // return x
|
|
|
513 res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
|
|
|
514 res.w[0] = y.w[0];
|
|
|
515 }
|
|
|
516 BID_RETURN (res)
|
|
|
517 } else { // at least one is infinity
|
|
|
518 if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
|
|
|
519 x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
|
|
|
520 x.w[0] = 0x0ull;
|
|
|
521 }
|
|
|
522 if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
|
|
|
523 y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
|
|
|
524 y.w[0] = 0x0ull;
|
|
|
525 }
|
|
|
526 }
|
|
|
527 }
|
|
|
528 // neither x nor y is NaN
|
|
|
529
|
|
|
530 // if not infinity, check for non-canonical values x (treated as zero)
|
|
|
531 if ((x.w[1] & MASK_ANY_INF) != MASK_INF) { // x != inf
|
|
|
532 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
|
|
|
533 // non-canonical
|
|
|
534 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
|
|
|
535 x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
|
|
|
536 x.w[0] = 0x0ull;
|
|
|
537 } else { // G0_G1 != 11
|
|
|
538 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
|
|
|
539 if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
|
|
|
540 ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
|
|
|
541 && x.w[0] > 0x378d8e63ffffffffull)) {
|
|
|
542 // x is non-canonical if coefficient is larger than 10^34 -1
|
|
|
543 x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
|
|
|
544 x.w[0] = 0x0ull;
|
|
|
545 } else { // canonical
|
|
|
546 ;
|
|
|
547 }
|
|
|
548 }
|
|
|
549 }
|
|
|
550 // no need to check for non-canonical y
|
|
|
551
|
|
|
552 // neither x nor y is NaN
|
|
|
553 tmp_fpsf = *pfpsf; // save fpsf
|
|
|
554 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
555 bid128_quiet_equal (&res1, &xnswp,
|
|
|
556 &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
557 _EXC_INFO_ARG);
|
|
|
558 bid128_quiet_greater (&res2, &xnswp,
|
|
|
559 &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
560 _EXC_INFO_ARG);
|
|
|
561 #else
|
|
|
562 res1 =
|
|
|
563 bid128_quiet_equal (xnswp,
|
|
|
564 ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
565 _EXC_INFO_ARG);
|
|
|
566 res2 =
|
|
|
567 bid128_quiet_greater (xnswp,
|
|
|
568 ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
569 _EXC_INFO_ARG);
|
|
|
570 #endif
|
|
|
571 *pfpsf = tmp_fpsf; // restore fpsf
|
|
|
572
|
|
|
573 if (res1) { // x = y
|
|
|
574 // return x with the sign of y
|
|
|
575 res.w[1] =
|
|
|
576 (x.w[1] & 0x7fffffffffffffffull) | (y.
|
|
|
577 w[1] & 0x8000000000000000ull);
|
|
|
578 res.w[0] = x.w[0];
|
|
|
579 } else if (res2) { // x > y
|
|
|
580 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
581 bid128_nextdown (&res,
|
|
|
582 &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
583 _EXC_INFO_ARG);
|
|
|
584 #else
|
|
|
585 res =
|
|
|
586 bid128_nextdown (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
587 _EXC_INFO_ARG);
|
|
|
588 #endif
|
|
|
589 BID_SWAP128 (res);
|
|
|
590 } else { // x < y
|
|
|
591 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
592 bid128_nextup (&res,
|
|
|
593 &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
|
594 #else
|
|
|
595 res =
|
|
|
596 bid128_nextup (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
|
597 #endif
|
|
|
598 BID_SWAP128 (res);
|
|
|
599 }
|
|
|
600 // if the operand x is finite but the result is infinite, signal
|
|
|
601 // overflow and inexact
|
|
|
602 if (((x.w[1] & MASK_SPECIAL) != MASK_SPECIAL)
|
|
|
603 && ((res.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
|
|
|
604 // set the inexact flag
|
|
|
605 *pfpsf |= INEXACT_EXCEPTION;
|
|
|
606 // set the overflow flag
|
|
|
607 *pfpsf |= OVERFLOW_EXCEPTION;
|
|
|
608 }
|
|
|
609 // if the result is in (-10^emin, 10^emin), and is different from the
|
|
|
610 // operand x, signal underflow and inexact
|
|
|
611 tmp1.w[HIGH_128W] = 0x0000314dc6448d93ull;
|
|
|
612 tmp1.w[LOW_128W] = 0x38c15b0a00000000ull; // +100...0[34] * 10^emin
|
|
|
613 tmp2.w[HIGH_128W] = res.w[1] & 0x7fffffffffffffffull;
|
|
|
614 tmp2.w[LOW_128W] = res.w[0];
|
|
|
615 tmp3.w[HIGH_128W] = res.w[1];
|
|
|
616 tmp3.w[LOW_128W] = res.w[0];
|
|
|
617 tmp_fpsf = *pfpsf; // save fpsf
|
|
|
618 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
619 bid128_quiet_greater (&res1, &tmp1,
|
|
|
620 &tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
621 _EXC_INFO_ARG);
|
|
|
622 bid128_quiet_not_equal (&res2, &xnswp,
|
|
|
623 &tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
624 _EXC_INFO_ARG);
|
|
|
625 #else
|
|
|
626 res1 =
|
|
|
627 bid128_quiet_greater (tmp1,
|
|
|
628 tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
629 _EXC_INFO_ARG);
|
|
|
630 res2 =
|
|
|
631 bid128_quiet_not_equal (xnswp,
|
|
|
632 tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
633 _EXC_INFO_ARG);
|
|
|
634 #endif
|
|
|
635 *pfpsf = tmp_fpsf; // restore fpsf
|
|
|
636 if (res1 && res2) {
|
|
|
637 // set the inexact flag
|
|
|
638 *pfpsf |= INEXACT_EXCEPTION;
|
|
|
639 // set the underflow flag
|
|
|
640 *pfpsf |= UNDERFLOW_EXCEPTION;
|
|
|
641 }
|
|
|
642 BID_RETURN (res);
|
|
|
643 }
|
