Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  1. /* Copyright (C) 2007-2015 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, 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_high (Qh, 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_high (Qh, 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, 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_high (Qh, 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_high (Qh, 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, 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_high (Qh, 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_high (Qh, 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, 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_high (Qh, 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_high (Qh, 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. }
  1852.