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