|
145
|
1 /* Copyright (C) 2007-2020 Free Software Foundation, Inc.
|
|
0
|
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 /*****************************************************************************
|
|
|
25 * BID64 add
|
|
|
26 *****************************************************************************
|
|
|
27 *
|
|
|
28 * Algorithm description:
|
|
|
29 *
|
|
|
30 * if(exponent_a < exponent_b)
|
|
|
31 * switch a, b
|
|
|
32 * diff_expon = exponent_a - exponent_b
|
|
|
33 * if(diff_expon > 16)
|
|
|
34 * return normalize(a)
|
|
|
35 * if(coefficient_a*10^diff_expon guaranteed below 2^62)
|
|
|
36 * S = sign_a*coefficient_a*10^diff_expon + sign_b*coefficient_b
|
|
|
37 * if(|S|<10^16)
|
|
|
38 * return get_BID64(sign(S),exponent_b,|S|)
|
|
|
39 * else
|
|
|
40 * determine number of extra digits in S (1, 2, or 3)
|
|
|
41 * return rounded result
|
|
|
42 * else // large exponent difference
|
|
|
43 * if(number_digits(coefficient_a*10^diff_expon) +/- 10^16)
|
|
|
44 * guaranteed the same as
|
|
|
45 * number_digits(coefficient_a*10^diff_expon) )
|
|
|
46 * S = normalize(coefficient_a + (sign_a^sign_b)*10^(16-diff_expon))
|
|
|
47 * corr = 10^16 + (sign_a^sign_b)*coefficient_b
|
|
|
48 * corr*10^exponent_b is rounded so it aligns with S*10^exponent_S
|
|
|
49 * return get_BID64(sign_a,exponent(S),S+rounded(corr))
|
|
|
50 * else
|
|
|
51 * add sign_a*coefficient_a*10^diff_expon, sign_b*coefficient_b
|
|
|
52 * in 128-bit integer arithmetic, then round to 16 decimal digits
|
|
|
53 *
|
|
|
54 *
|
|
|
55 ****************************************************************************/
|
|
|
56
|
|
|
57 #include "bid_internal.h"
|
|
|
58
|
|
|
59 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
60 void bid64_add (UINT64 * pres, UINT64 * px,
|
|
|
61 UINT64 *
|
|
|
62 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
|
|
63 _EXC_INFO_PARAM);
|
|
|
64 #else
|
|
|
65 UINT64 bid64_add (UINT64 x,
|
|
|
66 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
|
|
|
67 _EXC_MASKS_PARAM _EXC_INFO_PARAM);
|
|
|
68 #endif
|
|
|
69
|
|
|
70 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
71
|
|
|
72 void
|
|
|
73 bid64_sub (UINT64 * pres, UINT64 * px,
|
|
|
74 UINT64 *
|
|
|
75 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
|
|
76 _EXC_INFO_PARAM) {
|
|
|
77 UINT64 y = *py;
|
|
|
78 #if !DECIMAL_GLOBAL_ROUNDING
|
|
|
79 _IDEC_round rnd_mode = *prnd_mode;
|
|
|
80 #endif
|
|
|
81 // check if y is not NaN
|
|
|
82 if (((y & NAN_MASK64) != NAN_MASK64))
|
|
|
83 y ^= 0x8000000000000000ull;
|
|
|
84 bid64_add (pres, px,
|
|
|
85 &y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
86 _EXC_INFO_ARG);
|
|
|
87 }
|
|
|
88 #else
|
|
|
89
|
|
|
90 UINT64
|
|
|
91 bid64_sub (UINT64 x,
|
|
|
92 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
|
|
|
93 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
|
|
94 // check if y is not NaN
|
|
|
95 if (((y & NAN_MASK64) != NAN_MASK64))
|
|
|
96 y ^= 0x8000000000000000ull;
|
|
|
97
|
|
|
98 return bid64_add (x,
|
|
|
99 y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
|
100 _EXC_INFO_ARG);
|
|
|
101 }
|
|
|
102 #endif
|
|
|
103
|
|
|
104
|
|
|
105
|
|
|
106 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
107
|
|
|
108 void
|
|
|
109 bid64_add (UINT64 * pres, UINT64 * px,
|
|
|
110 UINT64 *
|
|
|
111 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
|
|
112 _EXC_INFO_PARAM) {
|
|
|
113 UINT64 x, y;
|
|
|
114 #else
|
|
|
115
|
|
|
116 UINT64
|
|
|
117 bid64_add (UINT64 x,
|
|
|
118 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
|
|
|
119 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
|
|
120 #endif
|
|
|
121
|
|
|
122 UINT128 CA, CT, CT_new;
|
|
|
123 UINT64 sign_x, sign_y, coefficient_x, coefficient_y, C64_new;
|
|
|
124 UINT64 valid_x, valid_y;
|
|
|
125 UINT64 res;
|
|
|
126 UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
|
|
|
127 rem_a;
|
|
|
128 UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp;
|
|
|
129 int_double tempx;
|
|
|
130 int exponent_x, exponent_y, exponent_a, exponent_b, diff_dec_expon;
|
|
|
131 int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
|
|
|
132 unsigned rmode, status;
|
|
|
133
|
|
|
134 #if DECIMAL_CALL_BY_REFERENCE
|
|
|
135 #if !DECIMAL_GLOBAL_ROUNDING
|
|
|
136 _IDEC_round rnd_mode = *prnd_mode;
|
|
|
137 #endif
|
|
|
138 x = *px;
|
|
|
139 y = *py;
|
|
|
140 #endif
|
|
|
141
|
|
|
142 valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
|
|
|
143 valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
|
|
|
144
|
|
|
145 // unpack arguments, check for NaN or Infinity
|
|
|
146 if (!valid_x) {
|
|
|
147 // x is Inf. or NaN
|
|
|
148
|
|
|
149 // test if x is NaN
|
|
|
150 if ((x & NAN_MASK64) == NAN_MASK64) {
|
|
|
151 #ifdef SET_STATUS_FLAGS
|
|
|
152 if (((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
|
|
|
153 || ((y & SNAN_MASK64) == SNAN_MASK64))
|
|
|
154 __set_status_flags (pfpsf, INVALID_EXCEPTION);
|
|
|
155 #endif
|
|
|
156 res = coefficient_x & QUIET_MASK64;
|
|
|
157 BID_RETURN (res);
|
|
|
158 }
|
|
|
159 // x is Infinity?
|
|
|
160 if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
|
|
|
161 // check if y is Inf
|
|
|
162 if (((y & NAN_MASK64) == INFINITY_MASK64)) {
|
|
|
163 if (sign_x == (y & 0x8000000000000000ull)) {
|
|
|
164 res = coefficient_x;
|
|
|
165 BID_RETURN (res);
|
|
|
166 }
|
|
|
167 // return NaN
|
|
|
168 {
|
|
|
169 #ifdef SET_STATUS_FLAGS
|
|
|
170 __set_status_flags (pfpsf, INVALID_EXCEPTION);
|
|
|
171 #endif
|
|
|
172 res = NAN_MASK64;
|
|
|
173 BID_RETURN (res);
|
|
|
174 }
|
|
|
175 }
|
|
|
176 // check if y is NaN
|
|
|
177 if (((y & NAN_MASK64) == NAN_MASK64)) {
|
|
|
178 res = coefficient_y & QUIET_MASK64;
|
|
|
179 #ifdef SET_STATUS_FLAGS
|
|
|
180 if (((y & SNAN_MASK64) == SNAN_MASK64))
|
|
|
181 __set_status_flags (pfpsf, INVALID_EXCEPTION);
|
|
|
182 #endif
|
|
|
183 BID_RETURN (res);
|
|
|
184 }
|
|
|
185 // otherwise return +/-Inf
|
|
|
186 {
|
|
|
187 res = coefficient_x;
|
|
|
188 BID_RETURN (res);
|
|
|
189 }
|
|
|
190 }
|
|
|
191 // x is 0
|
|
|
192 {
|
|
|
193 if (((y & INFINITY_MASK64) != INFINITY_MASK64) && coefficient_y) {
|
|
|
194 if (exponent_y <= exponent_x) {
|
|
|
195 res = y;
|
|
|
196 BID_RETURN (res);
|
|
|
197 }
|
|
|
198 }
|
|
|
199 }
|
|
|
200
|
|
|
201 }
|
|
|
202 if (!valid_y) {
|
|
|
203 // y is Inf. or NaN?
|
|
|
204 if (((y & INFINITY_MASK64) == INFINITY_MASK64)) {
|
|
|
205 #ifdef SET_STATUS_FLAGS
|
|
|
206 if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN
|
|
|
207 __set_status_flags (pfpsf, INVALID_EXCEPTION);
|
|
|
208 #endif
|
|
|
209 res = coefficient_y & QUIET_MASK64;
|
|
|
210 BID_RETURN (res);
|
|
|
211 }
|
|
|
212 // y is 0
|
|
|
213 if (!coefficient_x) { // x==0
|
|
|
214 if (exponent_x <= exponent_y)
|
|
|
215 res = ((UINT64) exponent_x) << 53;
|
|
|
216 else
|
|
|
217 res = ((UINT64) exponent_y) << 53;
|
|
|
218 if (sign_x == sign_y)
|
|
|
219 res |= sign_x;
|
|
|
220 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
|
|
221 #ifndef IEEE_ROUND_NEAREST
|
|
|
222 if (rnd_mode == ROUNDING_DOWN && sign_x != sign_y)
|
|
|
223 res |= 0x8000000000000000ull;
|
|
|
224 #endif
|
|
|
225 #endif
|
|
|
226 BID_RETURN (res);
|
|
|
227 } else if (exponent_y >= exponent_x) {
|
|
|
228 res = x;
|
|
|
229 BID_RETURN (res);
|
|
|
230 }
|
|
|
231 }
|
|
|
232 // sort arguments by exponent
|
|
|
233 if (exponent_x < exponent_y) {
|
|
|
234 sign_a = sign_y;
|
|
|
235 exponent_a = exponent_y;
|
|
|
236 coefficient_a = coefficient_y;
|
|
|
237 sign_b = sign_x;
|
|
|
238 exponent_b = exponent_x;
|
|
|
239 coefficient_b = coefficient_x;
|
|
|
240 } else {
|
|
|
241 sign_a = sign_x;
|
|
|
242 exponent_a = exponent_x;
|
|
|
243 coefficient_a = coefficient_x;
|
|
|
244 sign_b = sign_y;
|
|
|
245 exponent_b = exponent_y;
|
|
|
246 coefficient_b = coefficient_y;
|
|
|
247 }
|
|
|
248
|
|
|
249 // exponent difference
|
|
|
250 diff_dec_expon = exponent_a - exponent_b;
|
|
|
251
|
|
|
252 /* get binary coefficients of x and y */
|
|
|
253
|
|
|
254 //--- get number of bits in the coefficients of x and y ---
|
|
|
255
|
|
|
256 // version 2 (original)
|
|
|
257 tempx.d = (double) coefficient_a;
|
|
|
258 bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
|
|
|
259
|
|
|
260 if (diff_dec_expon > MAX_FORMAT_DIGITS) {
|
|
|
261 // normalize a to a 16-digit coefficient
|
|
|
262
|
|
|
263 scale_ca = estimate_decimal_digits[bin_expon_ca];
|
|
|
264 if (coefficient_a >= power10_table_128[scale_ca].w[0])
|
|
|
265 scale_ca++;
|
|
|
266
|
|
|
267 scale_k = 16 - scale_ca;
|
|
|
268
|
|
|
269 coefficient_a *= power10_table_128[scale_k].w[0];
|
|
|
270
|
|
|
271 diff_dec_expon -= scale_k;
|
|
|
272 exponent_a -= scale_k;
|
|
|
273
|
|
|
274 /* get binary coefficients of x and y */
|
|
|
275
|
|
|
276 //--- get number of bits in the coefficients of x and y ---
|
|
|
277 tempx.d = (double) coefficient_a;
|
|
|
278 bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
|
|
|
279
|
|
|
280 if (diff_dec_expon > MAX_FORMAT_DIGITS) {
|
|
|
281 #ifdef SET_STATUS_FLAGS
|
|
|
282 if (coefficient_b) {
|
|
|
283 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
|
|
284 }
|
|
|
285 #endif
|
|
|
286
|
|
|
287 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
|
|
288 #ifndef IEEE_ROUND_NEAREST
|
|
|
289 if (((rnd_mode) & 3) && coefficient_b) // not ROUNDING_TO_NEAREST
|
|
|
290 {
|
|
|
291 switch (rnd_mode) {
|
|
|
292 case ROUNDING_DOWN:
|
|
|
293 if (sign_b) {
|
|
|
294 coefficient_a -= ((((SINT64) sign_a) >> 63) | 1);
|
|
|
295 if (coefficient_a < 1000000000000000ull) {
|
|
|
296 exponent_a--;
|
|
|
297 coefficient_a = 9999999999999999ull;
|
|
|
298 } else if (coefficient_a >= 10000000000000000ull) {
|
|
|
299 exponent_a++;
|
|
|
300 coefficient_a = 1000000000000000ull;
|
|
|
301 }
|
|
|
302 }
|
|
|
303 break;
|
|
|
304 case ROUNDING_UP:
|
|
|
305 if (!sign_b) {
|
|
|
306 coefficient_a += ((((SINT64) sign_a) >> 63) | 1);
|
|
|
307 if (coefficient_a < 1000000000000000ull) {
|
|
|
308 exponent_a--;
|
|
|
309 coefficient_a = 9999999999999999ull;
|
|
|
310 } else if (coefficient_a >= 10000000000000000ull) {
|
|
|
311 exponent_a++;
|
|
|
312 coefficient_a = 1000000000000000ull;
|
|
|
313 }
|
|
|
314 }
|
|
|
315 break;
|
|
|
316 default: // RZ
|
|
|
317 if (sign_a != sign_b) {
|
|
|
318 coefficient_a--;
|
|
|
319 if (coefficient_a < 1000000000000000ull) {
|
|
|
320 exponent_a--;
|
|
|
321 coefficient_a = 9999999999999999ull;
|
|
|
322 }
|
|
|
323 }
|
|
|
324 break;
|
|
|
325 }
|
|
|
326 } else
|
|
|
327 #endif
|
|
|
328 #endif
|
|
|
329 // check special case here
|
|
|
330 if ((coefficient_a == 1000000000000000ull)
|
|
|
331 && (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
|
|
|
332 && (sign_a ^ sign_b)
|
|
|
333 && (coefficient_b > 5000000000000000ull)) {
|
|
|
334 coefficient_a = 9999999999999999ull;
|
|
|
335 exponent_a--;
|
|
|
336 }
|
|
|
337
|
|
|
338 res =
|
|
|
339 fast_get_BID64_check_OF (sign_a, exponent_a, coefficient_a,
|
|
|
340 rnd_mode, pfpsf);
|
|
|
341 BID_RETURN (res);
|
|
|
342 }
|
|
|
343 }
|
|
|
344 // test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62
|
|
|
345 if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) {
|
|
|
346 // coefficient_a*10^(exponent_a-exponent_b)<2^63
|
|
|
347
|
|
|
348 // multiply by 10^(exponent_a-exponent_b)
|
|
|
349 coefficient_a *= power10_table_128[diff_dec_expon].w[0];
|
|
|
350
|
|
|
351 // sign mask
|
|
|
352 sign_b = ((SINT64) sign_b) >> 63;
|
|
|
353 // apply sign to coeff. of b
|
|
|
354 coefficient_b = (coefficient_b + sign_b) ^ sign_b;
|
|
|
355
|
|
|
356 // apply sign to coefficient a
|
|
|
357 sign_a = ((SINT64) sign_a) >> 63;
|
|
|
358 coefficient_a = (coefficient_a + sign_a) ^ sign_a;
|
|
|
359
|
|
|
360 coefficient_a += coefficient_b;
|
|
|
361 // get sign
|
|
|
362 sign_s = ((SINT64) coefficient_a) >> 63;
|
|
|
363 coefficient_a = (coefficient_a + sign_s) ^ sign_s;
|
|
|
364 sign_s &= 0x8000000000000000ull;
|
|
|
365
|
|
|
366 // coefficient_a < 10^16 ?
|
|
|
367 if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) {
|
|
|
368 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
|
|
369 #ifndef IEEE_ROUND_NEAREST
|
|
|
370 if (rnd_mode == ROUNDING_DOWN && (!coefficient_a)
|
|
|
371 && sign_a != sign_b)
|
|
|
372 sign_s = 0x8000000000000000ull;
|
|
|
373 #endif
|
|
|
374 #endif
|
|
|
375 res = very_fast_get_BID64 (sign_s, exponent_b, coefficient_a);
|
|
|
376 BID_RETURN (res);
|
|
|
377 }
|
|
|
378 // otherwise rounding is necessary
|
|
|
379
|
|
|
380 // already know coefficient_a<10^19
|
|
|
381 // coefficient_a < 10^17 ?
|
|
|
382 if (coefficient_a < power10_table_128[17].w[0])
|
|
|
383 extra_digits = 1;
|
|
|
384 else if (coefficient_a < power10_table_128[18].w[0])
|
|
|
385 extra_digits = 2;
|
|
|
386 else
|
|
|
387 extra_digits = 3;
|
|
|
388
|
|
|
389 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
|
|
390 #ifndef IEEE_ROUND_NEAREST
|
|
|
391 rmode = rnd_mode;
|
|
|
392 if (sign_s && (unsigned) (rmode - 1) < 2)
|
|
|
393 rmode = 3 - rmode;
|
|
|
394 #else
|
|
|
395 rmode = 0;
|
|
|
396 #endif
|
|
|
397 #else
|
|
|
398 rmode = 0;
|
|
|
399 #endif
|
|
|
400 coefficient_a += round_const_table[rmode][extra_digits];
|
|
|
401
|
|
|
402 // get P*(2^M[extra_digits])/10^extra_digits
|
|
|
403 __mul_64x64_to_128 (CT, coefficient_a,
|
|
|
404 reciprocals10_64[extra_digits]);
|
|
|
405
|
|
|
406 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
|
|
407 amount = short_recip_scale[extra_digits];
|
|
|
408 C64 = CT.w[1] >> amount;
|
|
|
409
|
|
|
410 } else {
|
|
|
411 // coefficient_a*10^(exponent_a-exponent_b) is large
|
|
|
412 sign_s = sign_a;
|
|
|
413
|
|
|
414 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
|
|
415 #ifndef IEEE_ROUND_NEAREST
|
|
|
416 rmode = rnd_mode;
|
|
|
417 if (sign_s && (unsigned) (rmode - 1) < 2)
|
|
|
418 rmode = 3 - rmode;
|
|
|
419 #else
|
|
|
420 rmode = 0;
|
|
|
421 #endif
|
|
|
422 #else
|
|
|
423 rmode = 0;
|
|
|
424 #endif
|
|
|
425
|
|
|
426 // check whether we can take faster path
|
|
|
427 scale_ca = estimate_decimal_digits[bin_expon_ca];
|
|
|
428
|
|
|
429 sign_ab = sign_a ^ sign_b;
|
|
|
430 sign_ab = ((SINT64) sign_ab) >> 63;
|
|
|
431
|
|
|
432 // T1 = 10^(16-diff_dec_expon)
|
|
|
433 T1 = power10_table_128[16 - diff_dec_expon].w[0];
|
|
|
434
|
|
|
435 // get number of digits in coefficient_a
|
|
|
436 if (coefficient_a >= power10_table_128[scale_ca].w[0]) {
|
|
|
437 scale_ca++;
|
|
|
438 }
|
|
|
439
|
|
|
440 scale_k = 16 - scale_ca;
|
|
|
441
|
|
|
442 // addition
|
|
|
443 saved_ca = coefficient_a - T1;
|
|
|
444 coefficient_a =
|
|
|
445 (SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0];
|
|
|
446 extra_digits = diff_dec_expon - scale_k;
|
|
|
447
|
|
|
448 // apply sign
|
|
|
449 saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
|
|
|
450 // add 10^16 and rounding constant
|
|
|
451 coefficient_b =
|
|
|
452 saved_cb + 10000000000000000ull +
|
|
|
453 round_const_table[rmode][extra_digits];
|
|
|
454
|
|
|
455 // get P*(2^M[extra_digits])/10^extra_digits
|
|
|
456 __mul_64x64_to_128 (CT, coefficient_b,
|
|
|
457 reciprocals10_64[extra_digits]);
|
|
|
458
|
|
|
459 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
|
|
460 amount = short_recip_scale[extra_digits];
|
|
|
461 C0_64 = CT.w[1] >> amount;
|
|
|
462
|
|
|
463 // result coefficient
|
|
|
464 C64 = C0_64 + coefficient_a;
|
|
|
465 // filter out difficult (corner) cases
|
|
|
466 // this test ensures the number of digits in coefficient_a does not change
|
|
|
467 // after adding (the appropriately scaled and rounded) coefficient_b
|
|
|
468 if ((UINT64) (C64 - 1000000000000000ull - 1) >
|
|
|
469 9000000000000000ull - 2) {
|
|
|
470 if (C64 >= 10000000000000000ull) {
|
|
|
471 // result has more than 16 digits
|
|
|
472 if (!scale_k) {
|
|
|
473 // must divide coeff_a by 10
|
|
|
474 saved_ca = saved_ca + T1;
|
|
|
475 __mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull);
|
|
|
476 //reciprocals10_64[1]);
|
|
|
477 coefficient_a = CA.w[1] >> 1;
|
|
|
478 rem_a =
|
|
|
479 saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
|
|
|
480 coefficient_a = coefficient_a - T1;
|
|
|
481
|
|
|
482 saved_cb += rem_a * power10_table_128[diff_dec_expon].w[0];
|
|
|
483 } else
|
|
|
484 coefficient_a =
|
|
|
485 (SINT64) (saved_ca - T1 -
|
|
|
486 (T1 << 3)) * (SINT64) power10_table_128[scale_k -
|
|
|
487 1].w[0];
|
|
|
488
|
|
|
489 extra_digits++;
|
|
|
490 coefficient_b =
|
|
|
491 saved_cb + 100000000000000000ull +
|
|
|
492 round_const_table[rmode][extra_digits];
|
|
|
493
|
|
|
494 // get P*(2^M[extra_digits])/10^extra_digits
|
|
|
495 __mul_64x64_to_128 (CT, coefficient_b,
|
|
|
496 reciprocals10_64[extra_digits]);
|
|
|
497
|
|
|
498 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
|
|
499 amount = short_recip_scale[extra_digits];
|
|
|
500 C0_64 = CT.w[1] >> amount;
|
|
|
501
|
|
|
502 // result coefficient
|
|
|
503 C64 = C0_64 + coefficient_a;
|
|
|
504 } else if (C64 <= 1000000000000000ull) {
|
|
|
505 // less than 16 digits in result
|
|
|
506 coefficient_a =
|
|
|
507 (SINT64) saved_ca *(SINT64) power10_table_128[scale_k +
|
|
|
508 1].w[0];
|
|
|
509 //extra_digits --;
|
|
|
510 exponent_b--;
|
|
|
511 coefficient_b =
|
|
|
512 (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
|
|
|
513 round_const_table[rmode][extra_digits];
|
|
|
514
|
|
|
515 // get P*(2^M[extra_digits])/10^extra_digits
|
|
|
516 __mul_64x64_to_128 (CT_new, coefficient_b,
|
|
|
517 reciprocals10_64[extra_digits]);
|
|
|
518
|
|
|
519 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
|
|
520 amount = short_recip_scale[extra_digits];
|
|
|
521 C0_64 = CT_new.w[1] >> amount;
|
|
|
522
|
|
|
523 // result coefficient
|
|
|
524 C64_new = C0_64 + coefficient_a;
|
|
|
525 if (C64_new < 10000000000000000ull) {
|
|
|
526 C64 = C64_new;
|
|
|
527 #ifdef SET_STATUS_FLAGS
|
|
|
528 CT = CT_new;
|
|
|
529 #endif
|
|
|
530 } else
|
|
|
531 exponent_b++;
|
|
|
532 }
|
|
|
533
|
|
|
534 }
|
|
|
535
|
|
|
536 }
|
|
|
537
|
|
|
538 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
|
|
539 #ifndef IEEE_ROUND_NEAREST
|
|
|
540 if (rmode == 0) //ROUNDING_TO_NEAREST
|
|
|
541 #endif
|
|
|
542 if (C64 & 1) {
|
|
|
543 // check whether fractional part of initial_P/10^extra_digits is
|
|
|
544 // exactly .5
|
|
|
545 // this is the same as fractional part of
|
|
|
546 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
|
|
|
547
|
|
|
548 // get remainder
|
|
|
549 remainder_h = CT.w[1] << (64 - amount);
|
|
|
550
|
|
|
551 // test whether fractional part is 0
|
|
|
552 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
|
|
|
553 C64--;
|
|
|
554 }
|
|
|
555 }
|
|
|
556 #endif
|
|
|
557
|
|
|
558 #ifdef SET_STATUS_FLAGS
|
|
|
559 status = INEXACT_EXCEPTION;
|
|
|
560
|
|
|
561 // get remainder
|
|
|
562 remainder_h = CT.w[1] << (64 - amount);
|
|
|
563
|
|
|
564 switch (rmode) {
|
|
|
565 case ROUNDING_TO_NEAREST:
|
|
|
566 case ROUNDING_TIES_AWAY:
|
|
|
567 // test whether fractional part is 0
|
|
|
568 if ((remainder_h == 0x8000000000000000ull)
|
|
|
569 && (CT.w[0] < reciprocals10_64[extra_digits]))
|
|
|
570 status = EXACT_STATUS;
|
|
|
571 break;
|
|
|
572 case ROUNDING_DOWN:
|
|
|
573 case ROUNDING_TO_ZERO:
|
|
|
574 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
|
|
|
575 status = EXACT_STATUS;
|
|
|
576 //if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y;
|
|
|
577 break;
|
|
|
578 default:
|
|
|
579 // round up
|
|
|
580 __add_carry_out (tmp, carry, CT.w[0],
|
|
|
581 reciprocals10_64[extra_digits]);
|
|
|
582 if ((remainder_h >> (64 - amount)) + carry >=
|
|
|
583 (((UINT64) 1) << amount))
|
|
|
584 status = EXACT_STATUS;
|
|
|
585 break;
|
|
|
586 }
|
|
|
587 __set_status_flags (pfpsf, status);
|
|
|
588
|
|
|
589 #endif
|
|
|
590
|
|
|
591 res =
|
|
|
592 fast_get_BID64_check_OF (sign_s, exponent_b + extra_digits, C64,
|
|
|
593 rnd_mode, pfpsf);
|
|
|
594 BID_RETURN (res);
|
|
|
595 }
|
