Mercurial > hg > CbC > CbC_gcc

comparison libgcc/config/libbid/bid128_div.c @ 0:a06113de4d67

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first commit
author kent <kent@cr.ie.u-ryukyu.ac.jp>
date Fri, 17 Jul 2009 14:47:48 +0900
parents
children 04ced10e8804
comparison
equal deleted inserted replaced
-1:000000000000 0:a06113de4d67
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_div_macros.h"
26 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
27 #include <fenv.h>
28
29 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
30 #endif
31
32 extern UINT32 convert_table[5][128][2];
33 extern SINT8 factors[][2];
34 extern UINT8 packed_10000_zeros[];
35
36 BID128_FUNCTION_ARG2 (bid128_div, x, y)
37
38 UINT256 CA4, CA4r, P256;
39 UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
40 UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
41 valid_y;
42 int_float fx, fy, f64;
43 UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
44 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
45 digits_q, amount;
46 int nzeros, i, j, k, d5;
47 unsigned rmode;
48 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
49 fexcept_t binaryflags = 0;
50 #endif
51
52 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
53
54 // unpack arguments, check for NaN or Infinity
55 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
56 // test if x is NaN
57 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
58 #ifdef SET_STATUS_FLAGS
59 if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN
60 (y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull)
61 __set_status_flags (pfpsf, INVALID_EXCEPTION);
62 #endif
63 res.w[1] = (CX.w[1]) & QUIET_MASK64;
64 res.w[0] = CX.w[0];
65 BID_RETURN (res);
66 }
67 // x is Infinity?
68 if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
69 // check if y is Inf.
70 if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
71 // return NaN
72 {
73 #ifdef SET_STATUS_FLAGS
74 __set_status_flags (pfpsf, INVALID_EXCEPTION);
75 #endif
76 res.w[1] = 0x7c00000000000000ull;
77 res.w[0] = 0;
78 BID_RETURN (res);
79 }
80 // y is NaN?
81 if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull))
82 // return NaN
83 {
84 // return +/-Inf
85 res.w[1] = ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) |
86 0x7800000000000000ull;
87 res.w[0] = 0;
88 BID_RETURN (res);
89 }
90 }
91 // x is 0
92 if ((y.w[1] & 0x7800000000000000ull) < 0x7800000000000000ull) {
93 if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
94 #ifdef SET_STATUS_FLAGS
95 __set_status_flags (pfpsf, INVALID_EXCEPTION);
96 #endif
97 // x=y=0, return NaN
98 res.w[1] = 0x7c00000000000000ull;
99 res.w[0] = 0;
100 BID_RETURN (res);
101 }
102 // return 0
103 res.w[1] = (x.w[1] ^ y.w[1]) & 0x8000000000000000ull;
104 exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
105 if (exponent_x > DECIMAL_MAX_EXPON_128)
106 exponent_x = DECIMAL_MAX_EXPON_128;
107 else if (exponent_x < 0)
108 exponent_x = 0;
109 res.w[1] |= (((UINT64) exponent_x) << 49);
110 res.w[0] = 0;
111 BID_RETURN (res);
112 }
113 }
114 if (!valid_y) {
115 // y is Inf. or NaN
116
117 // test if y is NaN
118 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
119 #ifdef SET_STATUS_FLAGS
120 if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
121 __set_status_flags (pfpsf, INVALID_EXCEPTION);
122 #endif
123 res.w[1] = CY.w[1] & QUIET_MASK64;
124 res.w[0] = CY.w[0];
125 BID_RETURN (res);
126 }
127 // y is Infinity?
128 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
129 // return +/-0
130 res.w[1] = sign_x ^ sign_y;
131 res.w[0] = 0;
132 BID_RETURN (res);
133 }
134 // y is 0, return +/-Inf
135 #ifdef SET_STATUS_FLAGS
136 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
137 #endif
138 res.w[1] =
139 ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
140 res.w[0] = 0;
141 BID_RETURN (res);
142 }
143 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
144 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
145 #endif
146 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
147
148 if (__unsigned_compare_gt_128 (CY, CX)) {
149 // CX < CY
150
151 // 2^64
152 f64.i = 0x5f800000;
153
154 // fx ~ CX, fy ~ CY
155 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
156 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
157 // expon_cy - expon_cx
158 bin_index = (fy.i - fx.i) >> 23;
159
160 if (CX.w[1]) {
161 T = power10_index_binexp_128[bin_index].w[0];
162 __mul_64x128_short (CA, T, CX);
163 } else {
164 T128 = power10_index_binexp_128[bin_index];
165 __mul_64x128_short (CA, CX.w[0], T128);
166 }
167
168 ed2 = 33;
169 if (__unsigned_compare_gt_128 (CY, CA))
170 ed2++;
171
172 T128 = power10_table_128[ed2];
173 __mul_128x128_to_256 (CA4, CA, T128);
174
175 ed2 += estimate_decimal_digits[bin_index];
176 CQ.w[0] = CQ.w[1] = 0;
177 diff_expon = diff_expon - ed2;
178
179 } else {
180 // get CQ = CX/CY
181 __div_128_by_128 (&CQ, &CR, CX, CY);
182
183 if (!CR.w[1] && !CR.w[0]) {
184 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
185 pfpsf);
186 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
187 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
188 #endif
189 BID_RETURN (res);
190 }
191 // get number of decimal digits in CQ
192 // 2^64
193 f64.i = 0x5f800000;
194 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
195 // binary expon. of CQ
196 bin_expon = (fx.i - 0x3f800000) >> 23;
197
198 digits_q = estimate_decimal_digits[bin_expon];
199 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
200 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
201 if (__unsigned_compare_ge_128 (CQ, TP128))
202 digits_q++;
203
204 ed2 = 34 - digits_q;
205 T128.w[0] = power10_table_128[ed2].w[0];
206 T128.w[1] = power10_table_128[ed2].w[1];
207 __mul_128x128_to_256 (CA4, CR, T128);
208 diff_expon = diff_expon - ed2;
209 __mul_128x128_low (CQ, CQ, T128);
210
211 }
212
213 __div_256_by_128 (&CQ, &CA4, CY);
214
215 #ifdef SET_STATUS_FLAGS
216 if (CA4.w[0] || CA4.w[1]) {
217 // set status flags
218 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
219 }
220 #ifndef LEAVE_TRAILING_ZEROS
221 else
222 #endif
223 #else
224 #ifndef LEAVE_TRAILING_ZEROS
225 if (!CA4.w[0] && !CA4.w[1])
226 #endif
227 #endif
228 #ifndef LEAVE_TRAILING_ZEROS
229 // check whether result is exact
230 {
231 // check whether CX, CY are short
232 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
233 i = (int) CY.w[0] - 1;
234 j = (int) CX.w[0] - 1;
235 // difference in powers of 2 factors for Y and X
236 nzeros = ed2 - factors[i][0] + factors[j][0];
237 // difference in powers of 5 factors
238 d5 = ed2 - factors[i][1] + factors[j][1];
239 if (d5 < nzeros)
240 nzeros = d5;
241 // get P*(2^M[extra_digits])/10^extra_digits
242 __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
243
244 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
245 amount = recip_scale[nzeros];
246 __shr_128_long (CQ, Qh, amount);
247
248 diff_expon += nzeros;
249 } else {
250 // decompose Q as Qh*10^17 + Ql
251 //T128 = reciprocals10_128[17];
252 T128.w[0] = 0x44909befeb9fad49ull;
253 T128.w[1] = 0x000b877aa3236a4bull;
254 __mul_128x128_to_256 (P256, CQ, T128);
255 //amount = recip_scale[17];
256 Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
257 Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
258
259 if (!Q_low) {
260 diff_expon += 17;
261
262 tdigit[0] = Q_high & 0x3ffffff;
263 tdigit[1] = 0;
264 QX = Q_high >> 26;
265 QX32 = QX;
266 nzeros = 0;
267
268 for (j = 0; QX32; j++, QX32 >>= 7) {
269 k = (QX32 & 127);
270 tdigit[0] += convert_table[j][k][0];
271 tdigit[1] += convert_table[j][k][1];
272 if (tdigit[0] >= 100000000) {
273 tdigit[0] -= 100000000;
274 tdigit[1]++;
275 }
276 }
277
278 if (tdigit[1] >= 100000000) {
279 tdigit[1] -= 100000000;
280 if (tdigit[1] >= 100000000)
281 tdigit[1] -= 100000000;
282 }
283
284 digit = tdigit[0];
285 if (!digit && !tdigit[1])
286 nzeros += 16;
287 else {
288 if (!digit) {
289 nzeros += 8;
290 digit = tdigit[1];
291 }
292 // decompose digit
293 PD = (UINT64) digit *0x068DB8BBull;
294 digit_h = (UINT32) (PD >> 40);
295 digit_low = digit - digit_h * 10000;
296
297 if (!digit_low)
298 nzeros += 4;
299 else
300 digit_h = digit_low;
301
302 if (!(digit_h & 1))
303 nzeros +=
304 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
305 (digit_h & 7));
306 }
307
308 if (nzeros) {
309 __mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
310
311 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
312 amount = short_recip_scale[nzeros];
313 CQ.w[0] = CQ.w[1] >> amount;
314 } else
315 CQ.w[0] = Q_high;
316 CQ.w[1] = 0;
317
318 diff_expon += nzeros;
319 } else {
320 tdigit[0] = Q_low & 0x3ffffff;
321 tdigit[1] = 0;
322 QX = Q_low >> 26;
323 QX32 = QX;
324 nzeros = 0;
325
326 for (j = 0; QX32; j++, QX32 >>= 7) {
327 k = (QX32 & 127);
328 tdigit[0] += convert_table[j][k][0];
329 tdigit[1] += convert_table[j][k][1];
330 if (tdigit[0] >= 100000000) {
331 tdigit[0] -= 100000000;
332 tdigit[1]++;
333 }
334 }
335
336 if (tdigit[1] >= 100000000) {
337 tdigit[1] -= 100000000;
338 if (tdigit[1] >= 100000000)
339 tdigit[1] -= 100000000;
340 }
341
342 digit = tdigit[0];
343 if (!digit && !tdigit[1])
344 nzeros += 16;
345 else {
346 if (!digit) {
347 nzeros += 8;
348 digit = tdigit[1];
349 }
350 // decompose digit
351 PD = (UINT64) digit *0x068DB8BBull;
352 digit_h = (UINT32) (PD >> 40);
353 digit_low = digit - digit_h * 10000;
354
355 if (!digit_low)
356 nzeros += 4;
357 else
358 digit_h = digit_low;
359
360 if (!(digit_h & 1))
361 nzeros +=
362 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
363 (digit_h & 7));
364 }
365
366 if (nzeros) {
367 // get P*(2^M[extra_digits])/10^extra_digits
368 __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
369
370 //now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
371 amount = recip_scale[nzeros];
372 __shr_128 (CQ, Qh, amount);
373 }
374 diff_expon += nzeros;
375
376 }
377 }
378 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
379 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
380 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
381 #endif
382 BID_RETURN (res);
383 }
384 #endif
385
386 if (diff_expon >= 0) {
387 #ifdef IEEE_ROUND_NEAREST
388 // rounding
389 // 2*CA4 - CY
390 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
391 CA4r.w[0] = CA4.w[0] + CA4.w[0];
392 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
393 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
394
395 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
396 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
397
398 CQ.w[0] += carry64;
399 if (CQ.w[0] < carry64)
400 CQ.w[1]++;
401 #else
402 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
403 // rounding
404 // 2*CA4 - CY
405 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
406 CA4r.w[0] = CA4.w[0] + CA4.w[0];
407 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
408 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
409
410 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
411 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
412
413 CQ.w[0] += carry64;
414 if (CQ.w[0] < carry64)
415 CQ.w[1]++;
416 #else
417 rmode = rnd_mode;
418 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
419 rmode = 3 - rmode;
420 switch (rmode) {
421 case ROUNDING_TO_NEAREST: // round to nearest code
422 // rounding
423 // 2*CA4 - CY
424 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
425 CA4r.w[0] = CA4.w[0] + CA4.w[0];
426 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
427 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
428 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
429 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
430 CQ.w[0] += carry64;
431 if (CQ.w[0] < carry64)
432 CQ.w[1]++;
433 break;
434 case ROUNDING_TIES_AWAY:
435 // rounding
436 // 2*CA4 - CY
437 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
438 CA4r.w[0] = CA4.w[0] + CA4.w[0];
439 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
440 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
441 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
442 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
443 CQ.w[0] += carry64;
444 if (CQ.w[0] < carry64)
445 CQ.w[1]++;
446 break;
447 case ROUNDING_DOWN:
448 case ROUNDING_TO_ZERO:
449 break;
450 default: // rounding up
451 CQ.w[0]++;
452 if (!CQ.w[0])
453 CQ.w[1]++;
454 break;
455 }
456 #endif
457 #endif
458
459 } else {
460 #ifdef SET_STATUS_FLAGS
461 if (CA4.w[0] || CA4.w[1]) {
462 // set status flags
463 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
464 }
465 #endif
466
467 handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
468 CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
469 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
470 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
471 #endif
472 BID_RETURN (res);
473
474 }
475
476 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
477 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
478 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
479 #endif
480 BID_RETURN (res);
481 }
482
483
484 //#define LEAVE_TRAILING_ZEROS
485
486 TYPE0_FUNCTION_ARGTYPE1_ARGTYPE2 (UINT128, bid128dd_div, UINT64, x,
487 UINT64, y)
488
489 UINT256 CA4, CA4r, P256;
490 UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
491 UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
492 valid_y;
493 int_float fx, fy, f64;
494 UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
495 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
496 digits_q, amount;
497 int nzeros, i, j, k, d5;
498 unsigned rmode;
499 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
500 fexcept_t binaryflags = 0;
501 #endif
502
503 valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
504
505 // unpack arguments, check for NaN or Infinity
506 CX.w[1] = 0;
507 if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) {
508 #ifdef SET_STATUS_FLAGS
509 if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
510 __set_status_flags (pfpsf, INVALID_EXCEPTION);
511 #endif
512
513 // test if x is NaN
514 if ((x & NAN_MASK64) == NAN_MASK64) {
515 #ifdef SET_STATUS_FLAGS
516 if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
517 __set_status_flags (pfpsf, INVALID_EXCEPTION);
518 #endif
519 res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
520 __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
521 res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
522 BID_RETURN (res);
523 }
524 // x is Infinity?
525 if (((x) & 0x7800000000000000ull) == 0x7800000000000000ull) {
526 // check if y is Inf.
527 if ((((y) & 0x7c00000000000000ull) == 0x7800000000000000ull))
528 // return NaN
529 {
530 #ifdef SET_STATUS_FLAGS
531 __set_status_flags (pfpsf, INVALID_EXCEPTION);
532 #endif
533 res.w[1] = 0x7c00000000000000ull;
534 res.w[0] = 0;
535 BID_RETURN (res);
536 }
537 if ((((y) & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
538 // otherwise return +/-Inf
539 res.w[1] =
540 (((x) ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
541 res.w[0] = 0;
542 BID_RETURN (res);
543 }
544 }
545 // x is 0
546 if ((((y) & 0x7800000000000000ull) != 0x7800000000000000ull)) {
547 if(!CY.w[0]) {
548 #ifdef SET_STATUS_FLAGS
549 __set_status_flags (pfpsf, INVALID_EXCEPTION);
550 #endif
551 // x=y=0, return NaN
552 res.w[1] = 0x7c00000000000000ull;
553 res.w[0] = 0;
554 BID_RETURN (res);
555 }
556 // return 0
557 res.w[1] = ((x) ^ (y)) & 0x8000000000000000ull;
558 if (((y) & 0x6000000000000000ull) == 0x6000000000000000ull)
559 exponent_y = ((UINT32) ((y) >> 51)) & 0x3ff;
560 else
561 exponent_y = ((UINT32) ((y) >> 53)) & 0x3ff;
562 exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
563 if (exponent_x > DECIMAL_MAX_EXPON_128)
564 exponent_x = DECIMAL_MAX_EXPON_128;
565 else if (exponent_x < 0)
566 exponent_x = 0;
567 res.w[1] |= (((UINT64) exponent_x) << 49);
568 res.w[0] = 0;
569 BID_RETURN (res);
570 }
571 }
572
573 CY.w[1] = 0;
574 if (!valid_y) {
575 // y is Inf. or NaN
576
577 // test if y is NaN
578 if ((y & NAN_MASK64) == NAN_MASK64) {
579 #ifdef SET_STATUS_FLAGS
580 if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN
581 __set_status_flags (pfpsf, INVALID_EXCEPTION);
582 #endif
583 res.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
584 __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
585 res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
586 BID_RETURN (res);
587 }
588 // y is Infinity?
589 if (((y) & 0x7800000000000000ull) == 0x7800000000000000ull) {
590 // return +/-0
591 res.w[1] = sign_x ^ sign_y;
592 res.w[0] = 0;
593 BID_RETURN (res);
594 }
595 // y is 0, return +/-Inf
596 res.w[1] =
597 (((x) ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
598 res.w[0] = 0;
599 #ifdef SET_STATUS_FLAGS
600 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
601 #endif
602 BID_RETURN (res);
603 }
604 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
605 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
606 #endif
607 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
608
609 if (__unsigned_compare_gt_128 (CY, CX)) {
610 // CX < CY
611
612 // 2^64
613 f64.i = 0x5f800000;
614
615 // fx ~ CX, fy ~ CY
616 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
617 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
618 // expon_cy - expon_cx
619 bin_index = (fy.i - fx.i) >> 23;
620
621 if (CX.w[1]) {
622 T = power10_index_binexp_128[bin_index].w[0];
623 __mul_64x128_short (CA, T, CX);
624 } else {
625 T128 = power10_index_binexp_128[bin_index];
626 __mul_64x128_short (CA, CX.w[0], T128);
627 }
628
629 ed2 = 33;
630 if (__unsigned_compare_gt_128 (CY, CA))
631 ed2++;
632
633 T128 = power10_table_128[ed2];
634 __mul_128x128_to_256 (CA4, CA, T128);
635
636 ed2 += estimate_decimal_digits[bin_index];
637 CQ.w[0] = CQ.w[1] = 0;
638 diff_expon = diff_expon - ed2;
639
640 } else {
641 // get CQ = CX/CY
642 __div_128_by_128 (&CQ, &CR, CX, CY);
643
644 if (!CR.w[1] && !CR.w[0]) {
645 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
646 pfpsf);
647 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
648 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
649 #endif
650 BID_RETURN (res);
651 }
652 // get number of decimal digits in CQ
653 // 2^64
654 f64.i = 0x5f800000;
655 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
656 // binary expon. of CQ
657 bin_expon = (fx.i - 0x3f800000) >> 23;
658
659 digits_q = estimate_decimal_digits[bin_expon];
660 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
661 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
662 if (__unsigned_compare_ge_128 (CQ, TP128))
663 digits_q++;
664
665 ed2 = 34 - digits_q;
666 T128.w[0] = power10_table_128[ed2].w[0];
667 T128.w[1] = power10_table_128[ed2].w[1];
668 __mul_128x128_to_256 (CA4, CR, T128);
669 diff_expon = diff_expon - ed2;
670 __mul_128x128_low (CQ, CQ, T128);
671
672 }
673
674 __div_256_by_128 (&CQ, &CA4, CY);
675
676
677 #ifdef SET_STATUS_FLAGS
678 if (CA4.w[0] || CA4.w[1]) {
679 // set status flags
680 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
681 }
682 #ifndef LEAVE_TRAILING_ZEROS
683 else
684 #endif
685 #else
686 #ifndef LEAVE_TRAILING_ZEROS
687 if (!CA4.w[0] && !CA4.w[1])
688 #endif
689 #endif
690 #ifndef LEAVE_TRAILING_ZEROS
691 // check whether result is exact
692 {
693 // check whether CX, CY are short
694 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
695 i = (int) CY.w[0] - 1;
696 j = (int) CX.w[0] - 1;
697 // difference in powers of 2 factors for Y and X
698 nzeros = ed2 - factors[i][0] + factors[j][0];
699 // difference in powers of 5 factors
700 d5 = ed2 - factors[i][1] + factors[j][1];
701 if (d5 < nzeros)
702 nzeros = d5;
703 // get P*(2^M[extra_digits])/10^extra_digits
704 __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
705 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
706
707 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
708 amount = recip_scale[nzeros];
709 __shr_128_long (CQ, Qh, amount);
710
711 diff_expon += nzeros;
712 } else {
713 // decompose Q as Qh*10^17 + Ql
714 //T128 = reciprocals10_128[17];
715 T128.w[0] = 0x44909befeb9fad49ull;
716 T128.w[1] = 0x000b877aa3236a4bull;
717 __mul_128x128_to_256 (P256, CQ, T128);
718 //amount = recip_scale[17];
719 Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
720 Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
721
722 if (!Q_low) {
723 diff_expon += 17;
724
725 tdigit[0] = Q_high & 0x3ffffff;
726 tdigit[1] = 0;
727 QX = Q_high >> 26;
728 QX32 = QX;
729 nzeros = 0;
730
731 for (j = 0; QX32; j++, QX32 >>= 7) {
732 k = (QX32 & 127);
733 tdigit[0] += convert_table[j][k][0];
734 tdigit[1] += convert_table[j][k][1];
735 if (tdigit[0] >= 100000000) {
736 tdigit[0] -= 100000000;
737 tdigit[1]++;
738 }
739 }
740
741
742 if (tdigit[1] >= 100000000) {
743 tdigit[1] -= 100000000;
744 if (tdigit[1] >= 100000000)
745 tdigit[1] -= 100000000;
746 }
747
748 digit = tdigit[0];
749 if (!digit && !tdigit[1])
750 nzeros += 16;
751 else {
752 if (!digit) {
753 nzeros += 8;
754 digit = tdigit[1];
755 }
756 // decompose digit
757 PD = (UINT64) digit *0x068DB8BBull;
758 digit_h = (UINT32) (PD >> 40);
759 digit_low = digit - digit_h * 10000;
760
761 if (!digit_low)
762 nzeros += 4;
763 else
764 digit_h = digit_low;
765
766 if (!(digit_h & 1))
767 nzeros +=
768 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
769 (digit_h & 7));
770 }
771
772 if (nzeros) {
773 __mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
774
775 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
776 amount = short_recip_scale[nzeros];
777 CQ.w[0] = CQ.w[1] >> amount;
778 } else
779 CQ.w[0] = Q_high;
780 CQ.w[1] = 0;
781
782 diff_expon += nzeros;
783 } else {
784 tdigit[0] = Q_low & 0x3ffffff;
785 tdigit[1] = 0;
786 QX = Q_low >> 26;
787 QX32 = QX;
788 nzeros = 0;
789
790 for (j = 0; QX32; j++, QX32 >>= 7) {
791 k = (QX32 & 127);
792 tdigit[0] += convert_table[j][k][0];
793 tdigit[1] += convert_table[j][k][1];
794 if (tdigit[0] >= 100000000) {
795 tdigit[0] -= 100000000;
796 tdigit[1]++;
797 }
798 }
799
800 if (tdigit[1] >= 100000000) {
801 tdigit[1] -= 100000000;
802 if (tdigit[1] >= 100000000)
803 tdigit[1] -= 100000000;
804 }
805
806 digit = tdigit[0];
807 if (!digit && !tdigit[1])
808 nzeros += 16;
809 else {
810 if (!digit) {
811 nzeros += 8;
812 digit = tdigit[1];
813 }
814 // decompose digit
815 PD = (UINT64) digit *0x068DB8BBull;
816 digit_h = (UINT32) (PD >> 40);
817 digit_low = digit - digit_h * 10000;
818
819 if (!digit_low)
820 nzeros += 4;
821 else
822 digit_h = digit_low;
823
824 if (!(digit_h & 1))
825 nzeros +=
826 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
827 (digit_h & 7));
828 }
829
830 if (nzeros) {
831 // get P*(2^M[extra_digits])/10^extra_digits
832 __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
833
834 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
835 amount = recip_scale[nzeros];
836 __shr_128 (CQ, Qh, amount);
837 }
838 diff_expon += nzeros;
839
840 }
841 }
842 get_BID128(&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,pfpsf);
843 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
844 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
845 #endif
846 BID_RETURN (res);
847 }
848 #endif
849
850 if (diff_expon >= 0) {
851 #ifdef IEEE_ROUND_NEAREST
852 // rounding
853 // 2*CA4 - CY
854 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
855 CA4r.w[0] = CA4.w[0] + CA4.w[0];
856 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
857 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
858
859 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
860 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
861
862 CQ.w[0] += carry64;
863 if (CQ.w[0] < carry64)
864 CQ.w[1]++;
865 #else
866 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
867 // rounding
868 // 2*CA4 - CY
869 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
870 CA4r.w[0] = CA4.w[0] + CA4.w[0];
871 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
872 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
873
874 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
875 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
876
877 CQ.w[0] += carry64;
878 if (CQ.w[0] < carry64)
879 CQ.w[1]++;
880 #else
881 rmode = rnd_mode;
882 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
883 rmode = 3 - rmode;
884 switch (rmode) {
885 case ROUNDING_TO_NEAREST: // round to nearest code
886 // rounding
887 // 2*CA4 - CY
888 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
889 CA4r.w[0] = CA4.w[0] + CA4.w[0];
890 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
891 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
892 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
893 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
894 CQ.w[0] += carry64;
895 if (CQ.w[0] < carry64)
896 CQ.w[1]++;
897 break;
898 case ROUNDING_TIES_AWAY:
899 // rounding
900 // 2*CA4 - CY
901 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
902 CA4r.w[0] = CA4.w[0] + CA4.w[0];
903 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
904 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
905 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
906 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
907 CQ.w[0] += carry64;
908 if (CQ.w[0] < carry64)
909 CQ.w[1]++;
910 break;
911 case ROUNDING_DOWN:
912 case ROUNDING_TO_ZERO:
913 break;
914 default: // rounding up
915 CQ.w[0]++;
916 if (!CQ.w[0])
917 CQ.w[1]++;
918 break;
919 }
920 #endif
921 #endif
922
923 } else {
924 #ifdef SET_STATUS_FLAGS
925 if (CA4.w[0] || CA4.w[1]) {
926 // set status flags
927 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
928 }
929 #endif
930 handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
931 CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
932 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
933 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
934 #endif
935 BID_RETURN (res);
936
937 }
938
939 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
940 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
941 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
942 #endif
943 BID_RETURN (res);
944 }
945
946
947 BID128_FUNCTION_ARGTYPE1_ARG128 (bid128dq_div, UINT64, x, y)
948 UINT256 CA4, CA4r, P256;
949 UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
950 UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, valid_y,
951 PD;
952 int_float fx, fy, f64;
953 UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
954 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
955 digits_q, amount;
956 int nzeros, i, j, k, d5;
957 unsigned rmode;
958 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
959 fexcept_t binaryflags = 0;
960 #endif
961
962 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
963
964 // unpack arguments, check for NaN or Infinity
965 CX.w[1] = 0;
966 if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], x)) {
967 #ifdef SET_STATUS_FLAGS
968 if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
969 __set_status_flags (pfpsf, INVALID_EXCEPTION);
970 #endif
971
972 // test if x is NaN
973 if ((x & NAN_MASK64) == NAN_MASK64) {
974 #ifdef SET_STATUS_FLAGS
975 if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
976 __set_status_flags (pfpsf, INVALID_EXCEPTION);
977 #endif
978 res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
979 __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
980 res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
981 BID_RETURN (res);
982 }
983 // x is Infinity?
984 if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
985 // check if y is Inf.
986 if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
987 // return NaN
988 {
989 #ifdef SET_STATUS_FLAGS
990 __set_status_flags (pfpsf, INVALID_EXCEPTION);
991 #endif
992 res.w[1] = 0x7c00000000000000ull;
993 res.w[0] = 0;
994 BID_RETURN (res);
995 }
996 if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
997 // otherwise return +/-Inf
998 res.w[1] =
999 ((x ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
1000 res.w[0] = 0;
1001 BID_RETURN (res);
1002 }
1003 }
1004 // x is 0
1005 if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) {
1006 if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
1007 #ifdef SET_STATUS_FLAGS
1008 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1009 #endif
1010 // x=y=0, return NaN
1011 res.w[1] = 0x7c00000000000000ull;
1012 res.w[0] = 0;
1013 BID_RETURN (res);
1014 }
1015 // return 0
1016 res.w[1] = (x ^ y.w[1]) & 0x8000000000000000ull;
1017 exponent_x = exponent_x - exponent_y + (DECIMAL_EXPONENT_BIAS_128<<1) - DECIMAL_EXPONENT_BIAS;
1018 if (exponent_x > DECIMAL_MAX_EXPON_128)
1019 exponent_x = DECIMAL_MAX_EXPON_128;
1020 else if (exponent_x < 0)
1021 exponent_x = 0;
1022 res.w[1] |= (((UINT64) exponent_x) << 49);
1023 res.w[0] = 0;
1024 BID_RETURN (res);
1025 }
1026 }
1027 exponent_x += (DECIMAL_EXPONENT_BIAS_128 - DECIMAL_EXPONENT_BIAS);
1028
1029 if (!valid_y) {
1030 // y is Inf. or NaN
1031
1032 // test if y is NaN
1033 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
1034 #ifdef SET_STATUS_FLAGS
1035 if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
1036 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1037 #endif
1038 res.w[1] = CY.w[1] & QUIET_MASK64;
1039 res.w[0] = CY.w[0];
1040 BID_RETURN (res);
1041 }
1042 // y is Infinity?
1043 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
1044 // return +/-0
1045 res.w[1] = sign_x ^ sign_y;
1046 res.w[0] = 0;
1047 BID_RETURN (res);
1048 }
1049 // y is 0, return +/-Inf
1050 res.w[1] =
1051 ((x ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
1052 res.w[0] = 0;
1053 #ifdef SET_STATUS_FLAGS
1054 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
1055 #endif
1056 BID_RETURN (res);
1057 }
1058 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1059 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
1060 #endif
1061 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
1062
1063 if (__unsigned_compare_gt_128 (CY, CX)) {
1064 // CX < CY
1065
1066 // 2^64
1067 f64.i = 0x5f800000;
1068
1069 // fx ~ CX, fy ~ CY
1070 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
1071 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
1072 // expon_cy - expon_cx
1073 bin_index = (fy.i - fx.i) >> 23;
1074
1075 if (CX.w[1]) {
1076 T = power10_index_binexp_128[bin_index].w[0];
1077 __mul_64x128_short (CA, T, CX);
1078 } else {
1079 T128 = power10_index_binexp_128[bin_index];
1080 __mul_64x128_short (CA, CX.w[0], T128);
1081 }
1082
1083 ed2 = 33;
1084 if (__unsigned_compare_gt_128 (CY, CA))
1085 ed2++;
1086
1087 T128 = power10_table_128[ed2];
1088 __mul_128x128_to_256 (CA4, CA, T128);
1089
1090 ed2 += estimate_decimal_digits[bin_index];
1091 CQ.w[0] = CQ.w[1] = 0;
1092 diff_expon = diff_expon - ed2;
1093
1094 } else {
1095 // get CQ = CX/CY
1096 __div_128_by_128 (&CQ, &CR, CX, CY);
1097
1098 if (!CR.w[1] && !CR.w[0]) {
1099 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
1100 pfpsf);
1101 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1102 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1103 #endif
1104 BID_RETURN (res);
1105 }
1106 // get number of decimal digits in CQ
1107 // 2^64
1108 f64.i = 0x5f800000;
1109 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
1110 // binary expon. of CQ
1111 bin_expon = (fx.i - 0x3f800000) >> 23;
1112
1113 digits_q = estimate_decimal_digits[bin_expon];
1114 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
1115 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
1116 if (__unsigned_compare_ge_128 (CQ, TP128))
1117 digits_q++;
1118
1119 ed2 = 34 - digits_q;
1120 T128.w[0] = power10_table_128[ed2].w[0];
1121 T128.w[1] = power10_table_128[ed2].w[1];
1122 __mul_128x128_to_256 (CA4, CR, T128);
1123 diff_expon = diff_expon - ed2;
1124 __mul_128x128_low (CQ, CQ, T128);
1125
1126 }
1127
1128 __div_256_by_128 (&CQ, &CA4, CY);
1129
1130 #ifdef SET_STATUS_FLAGS
1131 if (CA4.w[0] || CA4.w[1]) {
1132 // set status flags
1133 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1134 }
1135 #ifndef LEAVE_TRAILING_ZEROS
1136 else
1137 #endif
1138 #else
1139 #ifndef LEAVE_TRAILING_ZEROS
1140 if (!CA4.w[0] && !CA4.w[1])
1141 #endif
1142 #endif
1143 #ifndef LEAVE_TRAILING_ZEROS
1144 // check whether result is exact
1145 {
1146 //printf("ed2=%d,nz=%d,a=%d,CQ="LX16","LX16", RH="LX16", RL="LX16"\n",ed2,nzeros,amount,CQ.w[1],CQ.w[0],reciprocals10_128[nzeros].w[1],reciprocals10_128[nzeros].w[0]);fflush(stdout);
1147 // check whether CX, CY are short
1148 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
1149 i = (int) CY.w[0] - 1;
1150 j = (int) CX.w[0] - 1;
1151 // difference in powers of 2 factors for Y and X
1152 nzeros = ed2 - factors[i][0] + factors[j][0];
1153 // difference in powers of 5 factors
1154 d5 = ed2 - factors[i][1] + factors[j][1];
1155 if (d5 < nzeros)
1156 nzeros = d5;
1157 // get P*(2^M[extra_digits])/10^extra_digits
1158 __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
1159 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1160
1161 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1162 amount = recip_scale[nzeros];
1163 __shr_128_long (CQ, Qh, amount);
1164
1165 diff_expon += nzeros;
1166 } else {
1167 // decompose Q as Qh*10^17 + Ql
1168 //T128 = reciprocals10_128[17];
1169 T128.w[0] = 0x44909befeb9fad49ull;
1170 T128.w[1] = 0x000b877aa3236a4bull;
1171 __mul_128x128_to_256 (P256, CQ, T128);
1172 //amount = recip_scale[17];
1173 Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
1174 Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
1175
1176 if (!Q_low) {
1177 diff_expon += 17;
1178
1179 tdigit[0] = Q_high & 0x3ffffff;
1180 tdigit[1] = 0;
1181 QX = Q_high >> 26;
1182 QX32 = QX;
1183 nzeros = 0;
1184
1185 for (j = 0; QX32; j++, QX32 >>= 7) {
1186 k = (QX32 & 127);
1187 tdigit[0] += convert_table[j][k][0];
1188 tdigit[1] += convert_table[j][k][1];
1189 if (tdigit[0] >= 100000000) {
1190 tdigit[0] -= 100000000;
1191 tdigit[1]++;
1192 }
1193 }
1194
1195
1196 if (tdigit[1] >= 100000000) {
1197 tdigit[1] -= 100000000;
1198 if (tdigit[1] >= 100000000)
1199 tdigit[1] -= 100000000;
1200 }
1201
1202 digit = tdigit[0];
1203 if (!digit && !tdigit[1])
1204 nzeros += 16;
1205 else {
1206 if (!digit) {
1207 nzeros += 8;
1208 digit = tdigit[1];
1209 }
1210 // decompose digit
1211 PD = (UINT64) digit *0x068DB8BBull;
1212 digit_h = (UINT32) (PD >> 40);
1213 //printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
1214 digit_low = digit - digit_h * 10000;
1215
1216 if (!digit_low)
1217 nzeros += 4;
1218 else
1219 digit_h = digit_low;
1220
1221 if (!(digit_h & 1))
1222 nzeros +=
1223 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1224 (digit_h & 7));
1225 }
1226
1227 if (nzeros) {
1228 __mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
1229
1230 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
1231 amount = short_recip_scale[nzeros];
1232 CQ.w[0] = CQ.w[1] >> amount;
1233 } else
1234 CQ.w[0] = Q_high;
1235 CQ.w[1] = 0;
1236
1237 diff_expon += nzeros;
1238 } else {
1239 tdigit[0] = Q_low & 0x3ffffff;
1240 tdigit[1] = 0;
1241 QX = Q_low >> 26;
1242 QX32 = QX;
1243 nzeros = 0;
1244
1245 for (j = 0; QX32; j++, QX32 >>= 7) {
1246 k = (QX32 & 127);
1247 tdigit[0] += convert_table[j][k][0];
1248 tdigit[1] += convert_table[j][k][1];
1249 if (tdigit[0] >= 100000000) {
1250 tdigit[0] -= 100000000;
1251 tdigit[1]++;
1252 }
1253 }
1254
1255 if (tdigit[1] >= 100000000) {
1256 tdigit[1] -= 100000000;
1257 if (tdigit[1] >= 100000000)
1258 tdigit[1] -= 100000000;
1259 }
1260
1261 digit = tdigit[0];
1262 if (!digit && !tdigit[1])
1263 nzeros += 16;
1264 else {
1265 if (!digit) {
1266 nzeros += 8;
1267 digit = tdigit[1];
1268 }
1269 // decompose digit
1270 PD = (UINT64) digit *0x068DB8BBull;
1271 digit_h = (UINT32) (PD >> 40);
1272 //printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
1273 digit_low = digit - digit_h * 10000;
1274
1275 if (!digit_low)
1276 nzeros += 4;
1277 else
1278 digit_h = digit_low;
1279
1280 if (!(digit_h & 1))
1281 nzeros +=
1282 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1283 (digit_h & 7));
1284 }
1285
1286 if (nzeros) {
1287 // get P*(2^M[extra_digits])/10^extra_digits
1288 __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
1289
1290 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1291 amount = recip_scale[nzeros];
1292 __shr_128 (CQ, Qh, amount);
1293 }
1294 diff_expon += nzeros;
1295
1296 }
1297 }
1298 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
1299 pfpsf);
1300 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1301 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1302 #endif
1303 BID_RETURN (res);
1304 }
1305 #endif
1306
1307 if (diff_expon >= 0) {
1308 #ifdef IEEE_ROUND_NEAREST
1309 // rounding
1310 // 2*CA4 - CY
1311 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1312 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1313 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1314 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1315
1316 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1317 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1318
1319 CQ.w[0] += carry64;
1320 if (CQ.w[0] < carry64)
1321 CQ.w[1]++;
1322 #else
1323 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1324 // rounding
1325 // 2*CA4 - CY
1326 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1327 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1328 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1329 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1330
1331 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1332 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1333
1334 CQ.w[0] += carry64;
1335 if (CQ.w[0] < carry64)
1336 CQ.w[1]++;
1337 #else
1338 rmode = rnd_mode;
1339 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
1340 rmode = 3 - rmode;
1341 switch (rmode) {
1342 case ROUNDING_TO_NEAREST: // round to nearest code
1343 // rounding
1344 // 2*CA4 - CY
1345 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1346 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1347 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1348 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1349 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1350 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1351 CQ.w[0] += carry64;
1352 if (CQ.w[0] < carry64)
1353 CQ.w[1]++;
1354 break;
1355 case ROUNDING_TIES_AWAY:
1356 // rounding
1357 // 2*CA4 - CY
1358 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1359 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1360 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1361 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1362 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1363 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1364 CQ.w[0] += carry64;
1365 if (CQ.w[0] < carry64)
1366 CQ.w[1]++;
1367 break;
1368 case ROUNDING_DOWN:
1369 case ROUNDING_TO_ZERO:
1370 break;
1371 default: // rounding up
1372 CQ.w[0]++;
1373 if (!CQ.w[0])
1374 CQ.w[1]++;
1375 break;
1376 }
1377 #endif
1378 #endif
1379
1380 } else {
1381 #ifdef SET_STATUS_FLAGS
1382 if (CA4.w[0] || CA4.w[1]) {
1383 // set status flags
1384 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1385 }
1386 #endif
1387 handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
1388 CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
1389 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1390 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1391 #endif
1392 BID_RETURN (res);
1393 }
1394
1395 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
1396 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1397 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1398 #endif
1399 BID_RETURN (res);
1400
1401 }
1402
1403
1404 BID128_FUNCTION_ARG128_ARGTYPE2 (bid128qd_div, x, UINT64, y)
1405 UINT256 CA4, CA4r, P256;
1406 UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
1407 UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
1408 valid_y;
1409 int_float fx, fy, f64;
1410 UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
1411 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
1412 digits_q, amount;
1413 int nzeros, i, j, k, d5, rmode;
1414 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1415 fexcept_t binaryflags = 0;
1416 #endif
1417
1418
1419 valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
1420 // unpack arguments, check for NaN or Infinity
1421 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
1422 // test if x is NaN
1423 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
1424 #ifdef SET_STATUS_FLAGS
1425 if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN
1426 (y & 0x7e00000000000000ull) == 0x7e00000000000000ull)
1427 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1428 #endif
1429 res.w[1] = (CX.w[1]) & QUIET_MASK64;
1430 res.w[0] = CX.w[0];
1431 BID_RETURN (res);
1432 }
1433 // x is Infinity?
1434 if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
1435 // check if y is Inf.
1436 if (((y & 0x7c00000000000000ull) == 0x7800000000000000ull))
1437 // return NaN
1438 {
1439 #ifdef SET_STATUS_FLAGS
1440 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1441 #endif
1442 res.w[1] = 0x7c00000000000000ull;
1443 res.w[0] = 0;
1444 BID_RETURN (res);
1445 }
1446 // y is NaN?
1447 if (((y & 0x7c00000000000000ull) != 0x7c00000000000000ull))
1448 // return NaN
1449 {
1450 // return +/-Inf
1451 res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull) |
1452 0x7800000000000000ull;
1453 res.w[0] = 0;
1454 BID_RETURN (res);
1455 }
1456 }
1457 // x is 0
1458 if ((y & 0x7800000000000000ull) < 0x7800000000000000ull) {
1459 if (!CY.w[0]) {
1460 #ifdef SET_STATUS_FLAGS
1461 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1462 #endif
1463 // x=y=0, return NaN
1464 res.w[1] = 0x7c00000000000000ull;
1465 res.w[0] = 0;
1466 BID_RETURN (res);
1467 }
1468 // return 0
1469 res.w[1] = (x.w[1] ^ y) & 0x8000000000000000ull;
1470 exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
1471 if (exponent_x > DECIMAL_MAX_EXPON_128)
1472 exponent_x = DECIMAL_MAX_EXPON_128;
1473 else if (exponent_x < 0)
1474 exponent_x = 0;
1475 res.w[1] |= (((UINT64) exponent_x) << 49);
1476 res.w[0] = 0;
1477 BID_RETURN (res);
1478 }
1479 }
1480 CY.w[1] = 0;
1481 if (!valid_y) {
1482 // y is Inf. or NaN
1483
1484 // test if y is NaN
1485 if ((y & NAN_MASK64) == NAN_MASK64) {
1486 #ifdef SET_STATUS_FLAGS
1487 if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN
1488 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1489 #endif
1490 res.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
1491 __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
1492 res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
1493 BID_RETURN (res);
1494 }
1495 // y is Infinity?
1496 if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
1497 // return +/-0
1498 res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull);
1499 res.w[0] = 0;
1500 BID_RETURN (res);
1501 }
1502 // y is 0
1503 #ifdef SET_STATUS_FLAGS
1504 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
1505 #endif
1506 res.w[1] = (sign_x ^ sign_y) | INFINITY_MASK64;
1507 res.w[0] = 0;
1508 BID_RETURN (res);
1509 }
1510 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1511 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
1512 #endif
1513 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
1514
1515 if (__unsigned_compare_gt_128 (CY, CX)) {
1516 // CX < CY
1517
1518 // 2^64
1519 f64.i = 0x5f800000;
1520
1521 // fx ~ CX, fy ~ CY
1522 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
1523 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
1524 // expon_cy - expon_cx
1525 bin_index = (fy.i - fx.i) >> 23;
1526
1527 if (CX.w[1]) {
1528 T = power10_index_binexp_128[bin_index].w[0];
1529 __mul_64x128_short (CA, T, CX);
1530 } else {
1531 T128 = power10_index_binexp_128[bin_index];
1532 __mul_64x128_short (CA, CX.w[0], T128);
1533 }
1534
1535 ed2 = 33;
1536 if (__unsigned_compare_gt_128 (CY, CA))
1537 ed2++;
1538
1539 T128 = power10_table_128[ed2];
1540 __mul_128x128_to_256 (CA4, CA, T128);
1541
1542 ed2 += estimate_decimal_digits[bin_index];
1543 CQ.w[0] = CQ.w[1] = 0;
1544 diff_expon = diff_expon - ed2;
1545
1546 } else {
1547 // get CQ = CX/CY
1548 __div_128_by_128 (&CQ, &CR, CX, CY);
1549
1550 if (!CR.w[1] && !CR.w[0]) {
1551 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
1552 pfpsf);
1553 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1554 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1555 #endif
1556 BID_RETURN (res);
1557 }
1558 // get number of decimal digits in CQ
1559 // 2^64
1560 f64.i = 0x5f800000;
1561 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
1562 // binary expon. of CQ
1563 bin_expon = (fx.i - 0x3f800000) >> 23;
1564
1565 digits_q = estimate_decimal_digits[bin_expon];
1566 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
1567 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
1568 if (__unsigned_compare_ge_128 (CQ, TP128))
1569 digits_q++;
1570
1571 ed2 = 34 - digits_q;
1572 T128.w[0] = power10_table_128[ed2].w[0];
1573 T128.w[1] = power10_table_128[ed2].w[1];
1574 __mul_128x128_to_256 (CA4, CR, T128);
1575 diff_expon = diff_expon - ed2;
1576 __mul_128x128_low (CQ, CQ, T128);
1577
1578 }
1579
1580 __div_256_by_128 (&CQ, &CA4, CY);
1581
1582
1583 #ifdef SET_STATUS_FLAGS
1584 if (CA4.w[0] || CA4.w[1]) {
1585 // set status flags
1586 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1587 }
1588 #ifndef LEAVE_TRAILING_ZEROS
1589 else
1590 #endif
1591 #else
1592 #ifndef LEAVE_TRAILING_ZEROS
1593 if (!CA4.w[0] && !CA4.w[1])
1594 #endif
1595 #endif
1596 #ifndef LEAVE_TRAILING_ZEROS
1597 // check whether result is exact
1598 {
1599 // check whether CX, CY are short
1600 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
1601 i = (int) CY.w[0] - 1;
1602 j = (int) CX.w[0] - 1;
1603 // difference in powers of 2 factors for Y and X
1604 nzeros = ed2 - factors[i][0] + factors[j][0];
1605 // difference in powers of 5 factors
1606 d5 = ed2 - factors[i][1] + factors[j][1];
1607 if (d5 < nzeros)
1608 nzeros = d5;
1609 // get P*(2^M[extra_digits])/10^extra_digits
1610 __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
1611 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1612
1613 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1614 amount = recip_scale[nzeros];
1615 __shr_128_long (CQ, Qh, amount);
1616
1617 diff_expon += nzeros;
1618 } else {
1619 // decompose Q as Qh*10^17 + Ql
1620 //T128 = reciprocals10_128[17];
1621 T128.w[0] = 0x44909befeb9fad49ull;
1622 T128.w[1] = 0x000b877aa3236a4bull;
1623 __mul_128x128_to_256 (P256, CQ, T128);
1624 //amount = recip_scale[17];
1625 Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
1626 Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
1627
1628 if (!Q_low) {
1629 diff_expon += 17;
1630
1631 tdigit[0] = Q_high & 0x3ffffff;
1632 tdigit[1] = 0;
1633 QX = Q_high >> 26;
1634 QX32 = QX;
1635 nzeros = 0;
1636
1637 for (j = 0; QX32; j++, QX32 >>= 7) {
1638 k = (QX32 & 127);
1639 tdigit[0] += convert_table[j][k][0];
1640 tdigit[1] += convert_table[j][k][1];
1641 if (tdigit[0] >= 100000000) {
1642 tdigit[0] -= 100000000;
1643 tdigit[1]++;
1644 }
1645 }
1646
1647
1648 if (tdigit[1] >= 100000000) {
1649 tdigit[1] -= 100000000;
1650 if (tdigit[1] >= 100000000)
1651 tdigit[1] -= 100000000;
1652 }
1653
1654 digit = tdigit[0];
1655 if (!digit && !tdigit[1])
1656 nzeros += 16;
1657 else {
1658 if (!digit) {
1659 nzeros += 8;
1660 digit = tdigit[1];
1661 }
1662 // decompose digit
1663 PD = (UINT64) digit *0x068DB8BBull;
1664 digit_h = (UINT32) (PD >> 40);
1665 digit_low = digit - digit_h * 10000;
1666
1667 if (!digit_low)
1668 nzeros += 4;
1669 else
1670 digit_h = digit_low;
1671
1672 if (!(digit_h & 1))
1673 nzeros +=
1674 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1675 (digit_h & 7));
1676 }
1677
1678 if (nzeros) {
1679 __mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
1680
1681 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
1682 amount = short_recip_scale[nzeros];
1683 CQ.w[0] = CQ.w[1] >> amount;
1684 } else
1685 CQ.w[0] = Q_high;
1686 CQ.w[1] = 0;
1687
1688 diff_expon += nzeros;
1689 } else {
1690 tdigit[0] = Q_low & 0x3ffffff;
1691 tdigit[1] = 0;
1692 QX = Q_low >> 26;
1693 QX32 = QX;
1694 nzeros = 0;
1695
1696 for (j = 0; QX32; j++, QX32 >>= 7) {
1697 k = (QX32 & 127);
1698 tdigit[0] += convert_table[j][k][0];
1699 tdigit[1] += convert_table[j][k][1];
1700 if (tdigit[0] >= 100000000) {
1701 tdigit[0] -= 100000000;
1702 tdigit[1]++;
1703 }
1704 }
1705
1706 if (tdigit[1] >= 100000000) {
1707 tdigit[1] -= 100000000;
1708 if (tdigit[1] >= 100000000)
1709 tdigit[1] -= 100000000;
1710 }
1711
1712 digit = tdigit[0];
1713 if (!digit && !tdigit[1])
1714 nzeros += 16;
1715 else {
1716 if (!digit) {
1717 nzeros += 8;
1718 digit = tdigit[1];
1719 }
1720 // decompose digit
1721 PD = (UINT64) digit *0x068DB8BBull;
1722 digit_h = (UINT32) (PD >> 40);
1723 digit_low = digit - digit_h * 10000;
1724
1725 if (!digit_low)
1726 nzeros += 4;
1727 else
1728 digit_h = digit_low;
1729
1730 if (!(digit_h & 1))
1731 nzeros +=
1732 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1733 (digit_h & 7));
1734 }
1735
1736 if (nzeros) {
1737 // get P*(2^M[extra_digits])/10^extra_digits
1738 __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
1739
1740 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1741 amount = recip_scale[nzeros];
1742 __shr_128 (CQ, Qh, amount);
1743 }
1744 diff_expon += nzeros;
1745
1746 }
1747 }
1748 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,pfpsf);
1749 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1750 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1751 #endif
1752 BID_RETURN (res);
1753 }
1754 #endif
1755
1756 if (diff_expon >= 0) {
1757 #ifdef IEEE_ROUND_NEAREST
1758 // rounding
1759 // 2*CA4 - CY
1760 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1761 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1762 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1763 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1764
1765 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1766 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1767
1768 CQ.w[0] += carry64;
1769 if (CQ.w[0] < carry64)
1770 CQ.w[1]++;
1771 #else
1772 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1773 // rounding
1774 // 2*CA4 - CY
1775 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1776 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1777 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1778 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1779
1780 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1781 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1782
1783 CQ.w[0] += carry64;
1784 if (CQ.w[0] < carry64)
1785 CQ.w[1]++;
1786 #else
1787 rmode = rnd_mode;
1788 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
1789 rmode = 3 - rmode;
1790 switch (rmode) {
1791 case ROUNDING_TO_NEAREST: // round to nearest code
1792 // rounding
1793 // 2*CA4 - CY
1794 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1795 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1796 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1797 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1798 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1799 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1800 CQ.w[0] += carry64;
1801 if (CQ.w[0] < carry64)
1802 CQ.w[1]++;
1803 break;
1804 case ROUNDING_TIES_AWAY:
1805 // rounding
1806 // 2*CA4 - CY
1807 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1808 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1809 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1810 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1811 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1812 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1813 CQ.w[0] += carry64;
1814 if (CQ.w[0] < carry64)
1815 CQ.w[1]++;
1816 break;
1817 case ROUNDING_DOWN:
1818 case ROUNDING_TO_ZERO:
1819 break;
1820 default: // rounding up
1821 CQ.w[0]++;
1822 if (!CQ.w[0])
1823 CQ.w[1]++;
1824 break;
1825 }
1826 #endif
1827 #endif
1828
1829 } else {
1830 #ifdef SET_STATUS_FLAGS
1831 if (CA4.w[0] || CA4.w[1]) {
1832 // set status flags
1833 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1834 }
1835 #endif
1836 handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
1837 CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
1838 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1839 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1840 #endif
1841 BID_RETURN (res);
1842
1843 }
1844
1845 get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
1846 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1847 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1848 #endif
1849 BID_RETURN (res);
1850
1851 }