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