OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-stable/] [gcc-4.5.1/] [libgcc/] [config/] [libbid/] [bid128_div.c] - Blame information for rev 853

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 272 jeremybenn 
/* 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 
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.