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. #ifndef __BIDECIMAL_H
  25. #define __BIDECIMAL_H
  26.  
  27. #include "bid_conf.h"
  28. #include "bid_functions.h"
  29.  
  30. #define __BID_INLINE__ static __inline
  31.  
  32. /*********************************************************************
  33.  *
  34.  *      Logical Shift Macros
  35.  *
  36.  *********************************************************************/
  37.  
  38. #define __shr_128(Q, A, k)                        \
  39. {                                                 \
  40.      (Q).w[0] = (A).w[0] >> k;                      \
  41.          (Q).w[0] |= (A).w[1] << (64-k);               \
  42.          (Q).w[1] = (A).w[1] >> k;                    \
  43. }
  44.  
  45. #define __shr_128_long(Q, A, k)                   \
  46. {                                                 \
  47.         if((k)<64) {                                  \
  48.      (Q).w[0] = (A).w[0] >> k;                    \
  49.          (Q).w[0] |= (A).w[1] << (64-k);              \
  50.          (Q).w[1] = (A).w[1] >> k;                    \
  51.         }                                             \
  52.         else {                                        \
  53.          (Q).w[0] = (A).w[1]>>((k)-64);               \
  54.          (Q).w[1] = 0;                                \
  55.         }                                             \
  56. }
  57.  
  58. #define __shl_128_long(Q, A, k)                   \
  59. {                                                 \
  60.         if((k)<64) {                                  \
  61.      (Q).w[1] = (A).w[1] << k;                    \
  62.          (Q).w[1] |= (A).w[0] >> (64-k);              \
  63.          (Q).w[0] = (A).w[0] << k;                    \
  64.         }                                             \
  65.         else {                                        \
  66.          (Q).w[1] = (A).w[0]<<((k)-64);               \
  67.          (Q).w[0] = 0;                                \
  68.         }                                             \
  69. }
  70.  
  71. #define __low_64(Q)  (Q).w[0]
  72. /*********************************************************************
  73.  *
  74.  *      String Macros
  75.  *
  76.  *********************************************************************/
  77. #define tolower_macro(x) (((unsigned char)((x)-'A')<=('Z'-'A'))?((x)-'A'+'a'):(x))
  78. /*********************************************************************
  79.  *
  80.  *      Compare Macros
  81.  *
  82.  *********************************************************************/
  83. // greater than
  84. //  return 0 if A<=B
  85. //  non-zero if A>B
  86. #define __unsigned_compare_gt_128(A, B)  \
  87.     ((A.w[1]>B.w[1]) || ((A.w[1]==B.w[1]) && (A.w[0]>B.w[0])))
  88. // greater-or-equal
  89. #define __unsigned_compare_ge_128(A, B)  \
  90.     ((A.w[1]>B.w[1]) || ((A.w[1]==B.w[1]) && (A.w[0]>=B.w[0])))
  91. #define __test_equal_128(A, B)  (((A).w[1]==(B).w[1]) && ((A).w[0]==(B).w[0]))
  92. // tighten exponent range
  93. #define __tight_bin_range_128(bp, P, bin_expon)  \
  94. {                                                \
  95. UINT64 M;                                        \
  96.         M = 1;                                       \
  97.         (bp) = (bin_expon);                          \
  98.         if((bp)<63) {                                \
  99.           M <<= ((bp)+1);                            \
  100.           if((P).w[0] >= M) (bp)++; }                 \
  101.         else if((bp)>64) {                           \
  102.           M <<= ((bp)+1-64);                         \
  103.           if(((P).w[1]>M) ||((P).w[1]==M && (P).w[0]))\
  104.               (bp)++; }                              \
  105.         else if((P).w[1]) (bp)++;                    \
  106. }
  107. /*********************************************************************
  108.  *
  109.  *      Add/Subtract Macros
  110.  *
  111.  *********************************************************************/
  112. // add 64-bit value to 128-bit
  113. #define __add_128_64(R128, A128, B64)    \
  114. {                                        \
  115. UINT64 R64H;                             \
  116.         R64H = (A128).w[1];                 \
  117.         (R128).w[0] = (B64) + (A128).w[0];     \
  118.         if((R128).w[0] < (B64))               \
  119.           R64H ++;                           \
  120.     (R128).w[1] = R64H;                  \
  121. }
  122. // subtract 64-bit value from 128-bit
  123. #define __sub_128_64(R128, A128, B64)    \
  124. {                                        \
  125. UINT64 R64H;                             \
  126.         R64H = (A128).w[1];                  \
  127.         if((A128).w[0] < (B64))               \
  128.           R64H --;                           \
  129.     (R128).w[1] = R64H;                  \
  130.         (R128).w[0] = (A128).w[0] - (B64);     \
  131. }
  132. // add 128-bit value to 128-bit
  133. // assume no carry-out
  134. #define __add_128_128(R128, A128, B128)  \
  135. {                                        \
  136. UINT128 Q128;                            \
  137.         Q128.w[1] = (A128).w[1]+(B128).w[1]; \
  138.         Q128.w[0] = (B128).w[0] + (A128).w[0];  \
  139.         if(Q128.w[0] < (B128).w[0])            \
  140.           Q128.w[1] ++;                      \
  141.     (R128).w[1] = Q128.w[1];             \
  142.     (R128).w[0] = Q128.w[0];               \
  143. }
  144. #define __sub_128_128(R128, A128, B128)  \
  145. {                                        \
  146. UINT128 Q128;                            \
  147.         Q128.w[1] = (A128).w[1]-(B128).w[1]; \
  148.         Q128.w[0] = (A128).w[0] - (B128).w[0];  \
  149.         if((A128).w[0] < (B128).w[0])          \
  150.           Q128.w[1] --;                      \
  151.     (R128).w[1] = Q128.w[1];             \
  152.     (R128).w[0] = Q128.w[0];               \
  153. }
  154. #define __add_carry_out(S, CY, X, Y)    \
  155. {                                      \
  156. UINT64 X1=X;                           \
  157.         S = X + Y;                         \
  158.         CY = (S<X1) ? 1 : 0;                \
  159. }
  160. #define __add_carry_in_out(S, CY, X, Y, CI)    \
  161. {                                             \
  162. UINT64 X1;                                    \
  163.         X1 = X + CI;                              \
  164.         S = X1 + Y;                               \
  165.         CY = ((S<X1) || (X1<CI)) ? 1 : 0;          \
  166. }
  167. #define __sub_borrow_out(S, CY, X, Y)    \
  168. {                                      \
  169. UINT64 X1=X;                           \
  170.         S = X - Y;                         \
  171.         CY = (S>X1) ? 1 : 0;                \
  172. }
  173. #define __sub_borrow_in_out(S, CY, X, Y, CI)    \
  174. {                                             \
  175. UINT64 X1, X0=X;                              \
  176.         X1 = X - CI;                              \
  177.         S = X1 - Y;                               \
  178.         CY = ((S>X1) || (X1>X0)) ? 1 : 0;          \
  179. }
  180. // increment C128 and check for rounding overflow:
  181. // if (C_128) = 10^34 then (C_128) = 10^33 and increment the exponent
  182. #define INCREMENT(C_128, exp)                                           \
  183. {                                                                       \
  184.   C_128.w[0]++;                                                         \
  185.   if (C_128.w[0] == 0) C_128.w[1]++;                                    \
  186.   if (C_128.w[1] == 0x0001ed09bead87c0ull &&                            \
  187.       C_128.w[0] == 0x378d8e6400000000ull) {                            \
  188.     exp++;                                                              \
  189.     C_128.w[1] = 0x0000314dc6448d93ull;                                 \
  190.     C_128.w[0] = 0x38c15b0a00000000ull;                                 \
  191.   }                                                                     \
  192. }
  193. // decrement C128 and check for rounding underflow, but only at the
  194. // boundary: if C_128 = 10^33 - 1 and exp > 0 then C_128 = 10^34 - 1
  195. // and decrement the exponent
  196. #define DECREMENT(C_128, exp)                                           \
  197. {                                                                       \
  198.   C_128.w[0]--;                                                         \
  199.   if (C_128.w[0] == 0xffffffffffffffffull) C_128.w[1]--;                \
  200.   if (C_128.w[1] == 0x0000314dc6448d93ull &&                            \
  201.       C_128.w[0] == 0x38c15b09ffffffffull && exp > 0) {                 \
  202.     exp--;                                                              \
  203.     C_128.w[1] = 0x0001ed09bead87c0ull;                                 \
  204.     C_128.w[0] = 0x378d8e63ffffffffull;                                 \
  205.   }                                                                     \
  206. }
  207.  
  208.  /*********************************************************************
  209.  *
  210.  *      Multiply Macros
  211.  *
  212.  *********************************************************************/
  213. #define __mul_64x64_to_64(P64, CX, CY)  (P64) = (CX) * (CY)
  214. /***************************************
  215.  *  Signed, Full 64x64-bit Multiply
  216.  ***************************************/
  217. #define __imul_64x64_to_128(P, CX, CY)  \
  218. {                                       \
  219. UINT64 SX, SY;                          \
  220.    __mul_64x64_to_128(P, CX, CY);       \
  221.                                         \
  222.    SX = ((SINT64)(CX))>>63;             \
  223.    SY = ((SINT64)(CY))>>63;             \
  224.    SX &= CY;   SY &= CX;                \
  225.                                         \
  226.    (P).w[1] = (P).w[1] - SX - SY;       \
  227. }
  228. /***************************************
  229.  *  Signed, Full 64x128-bit Multiply
  230.  ***************************************/
  231. #define __imul_64x128_full(Ph, Ql, A, B)          \
  232. {                                                 \
  233. UINT128 ALBL, ALBH, QM2, QM;                      \
  234.                                                   \
  235.         __imul_64x64_to_128(ALBH, (A), (B).w[1]);     \
  236.         __imul_64x64_to_128(ALBL, (A), (B).w[0]);      \
  237.                                                   \
  238.         (Ql).w[0] = ALBL.w[0];                          \
  239.         QM.w[0] = ALBL.w[1];                           \
  240.         QM.w[1] = ((SINT64)ALBL.w[1])>>63;            \
  241.     __add_128_128(QM2, ALBH, QM);                 \
  242.         (Ql).w[1] = QM2.w[0];                          \
  243.     Ph = QM2.w[1];                                \
  244. }
  245. /*****************************************************
  246.  *      Unsigned Multiply Macros
  247.  *****************************************************/
  248. // get full 64x64bit product
  249. //
  250. #define __mul_64x64_to_128(P, CX, CY)   \
  251. {                                       \
  252. UINT64 CXH, CXL, CYH,CYL,PL,PH,PM,PM2;\
  253.         CXH = (CX) >> 32;                     \
  254.         CXL = (UINT32)(CX);                   \
  255.         CYH = (CY) >> 32;                     \
  256.         CYL = (UINT32)(CY);                   \
  257.                                               \
  258.     PM = CXH*CYL;                         \
  259.         PH = CXH*CYH;                         \
  260.         PL = CXL*CYL;                         \
  261.         PM2 = CXL*CYH;                        \
  262.         PH += (PM>>32);                       \
  263.         PM = (UINT64)((UINT32)PM)+PM2+(PL>>32); \
  264.                                           \
  265.         (P).w[1] = PH + (PM>>32);             \
  266.         (P).w[0] = (PM<<32)+(UINT32)PL;       \
  267. }
  268. // get full 64x64bit product
  269. // Note:
  270. // This macro is used for CX < 2^61, CY < 2^61
  271. //
  272. #define __mul_64x64_to_128_fast(P, CX, CY)   \
  273. {                                       \
  274. UINT64 CXH, CXL, CYH, CYL, PL, PH, PM;  \
  275.         CXH = (CX) >> 32;                   \
  276.         CXL = (UINT32)(CX);                 \
  277.         CYH = (CY) >> 32;                   \
  278.         CYL = (UINT32)(CY);                 \
  279.                                             \
  280.     PM = CXH*CYL;                       \
  281.         PL = CXL*CYL;                       \
  282.         PH = CXH*CYH;                       \
  283.         PM += CXL*CYH;                      \
  284.         PM += (PL>>32);                     \
  285.                                         \
  286.         (P).w[1] = PH + (PM>>32);           \
  287.         (P).w[0] = (PM<<32)+(UINT32)PL;      \
  288. }
  289. // used for CX< 2^60
  290. #define __sqr64_fast(P, CX)   \
  291. {                                       \
  292. UINT64 CXH, CXL, PL, PH, PM;            \
  293.         CXH = (CX) >> 32;                   \
  294.         CXL = (UINT32)(CX);                 \
  295.                                             \
  296.     PM = CXH*CXL;                       \
  297.         PL = CXL*CXL;                       \
  298.         PH = CXH*CXH;                       \
  299.         PM += PM;                           \
  300.         PM += (PL>>32);                     \
  301.                                         \
  302.         (P).w[1] = PH + (PM>>32);           \
  303.         (P).w[0] = (PM<<32)+(UINT32)PL;     \
  304. }
  305. // get full 64x64bit product
  306. // Note:
  307. // This implementation is used for CX < 2^61, CY < 2^61
  308. //
  309. #define __mul_64x64_to_64_high_fast(P, CX, CY)   \
  310. {                                       \
  311. UINT64 CXH, CXL, CYH, CYL, PL, PH, PM;  \
  312.         CXH = (CX) >> 32;                   \
  313.         CXL = (UINT32)(CX);                 \
  314.         CYH = (CY) >> 32;                   \
  315.         CYL = (UINT32)(CY);                 \
  316.                                             \
  317.     PM = CXH*CYL;                       \
  318.         PL = CXL*CYL;                       \
  319.         PH = CXH*CYH;                       \
  320.         PM += CXL*CYH;                      \
  321.         PM += (PL>>32);                     \
  322.                                         \
  323.         (P) = PH + (PM>>32);                \
  324. }
  325. // get full 64x64bit product
  326. //
  327. #define __mul_64x64_to_128_full(P, CX, CY)     \
  328. {                                         \
  329. UINT64 CXH, CXL, CYH,CYL,PL,PH,PM,PM2;\
  330.         CXH = (CX) >> 32;                     \
  331.         CXL = (UINT32)(CX);                   \
  332.         CYH = (CY) >> 32;                     \
  333.         CYL = (UINT32)(CY);                   \
  334.                                               \
  335.     PM = CXH*CYL;                         \
  336.         PH = CXH*CYH;                         \
  337.         PL = CXL*CYL;                         \
  338.         PM2 = CXL*CYH;                        \
  339.         PH += (PM>>32);                       \
  340.         PM = (UINT64)((UINT32)PM)+PM2+(PL>>32); \
  341.                                           \
  342.         (P).w[1] = PH + (PM>>32);             \
  343.         (P).w[0] = (PM<<32)+(UINT32)PL;        \
  344. }
  345. #define __mul_128x128_high(Q, A, B)               \
  346. {                                                 \
  347. UINT128 ALBL, ALBH, AHBL, AHBH, QM, QM2;          \
  348.                                                   \
  349.         __mul_64x64_to_128(ALBH, (A).w[0], (B).w[1]);  \
  350.         __mul_64x64_to_128(AHBL, (B).w[0], (A).w[1]);  \
  351.         __mul_64x64_to_128(ALBL, (A).w[0], (B).w[0]);   \
  352.         __mul_64x64_to_128(AHBH, (A).w[1],(B).w[1]);  \
  353.                                                   \
  354.     __add_128_128(QM, ALBH, AHBL);                \
  355.     __add_128_64(QM2, QM, ALBL.w[1]);             \
  356.     __add_128_64((Q), AHBH, QM2.w[1]);            \
  357. }
  358. #define __mul_128x128_full(Qh, Ql, A, B)          \
  359. {                                                 \
  360. UINT128 ALBL, ALBH, AHBL, AHBH, QM, QM2;          \
  361.                                                   \
  362.         __mul_64x64_to_128(ALBH, (A).w[0], (B).w[1]);  \
  363.         __mul_64x64_to_128(AHBL, (B).w[0], (A).w[1]);  \
  364.         __mul_64x64_to_128(ALBL, (A).w[0], (B).w[0]);   \
  365.         __mul_64x64_to_128(AHBH, (A).w[1],(B).w[1]);  \
  366.                                                   \
  367.     __add_128_128(QM, ALBH, AHBL);                \
  368.         (Ql).w[0] = ALBL.w[0];                          \
  369.     __add_128_64(QM2, QM, ALBL.w[1]);             \
  370.     __add_128_64((Qh), AHBH, QM2.w[1]);           \
  371.         (Ql).w[1] = QM2.w[0];                          \
  372. }
  373. #define __mul_128x128_low(Ql, A, B)               \
  374. {                                                 \
  375. UINT128 ALBL;                                     \
  376. UINT64 QM64;                                      \
  377.                                                   \
  378.         __mul_64x64_to_128(ALBL, (A).w[0], (B).w[0]);   \
  379.         QM64 = (B).w[0]*(A).w[1] + (A).w[0]*(B).w[1];   \
  380.                                                   \
  381.         (Ql).w[0] = ALBL.w[0];                          \
  382.         (Ql).w[1] = QM64 + ALBL.w[1];                 \
  383. }
  384. #define __mul_64x128_low(Ql, A, B)                \
  385. {                                                 \
  386.   UINT128 ALBL, ALBH, QM2;                        \
  387.   __mul_64x64_to_128(ALBH, (A), (B).w[1]);        \
  388.   __mul_64x64_to_128(ALBL, (A), (B).w[0]);        \
  389.   (Ql).w[0] = ALBL.w[0];                          \
  390.   __add_128_64(QM2, ALBH, ALBL.w[1]);             \
  391.   (Ql).w[1] = QM2.w[0];                           \
  392. }
  393. #define __mul_64x128_full(Ph, Ql, A, B)           \
  394. {                                                 \
  395. UINT128 ALBL, ALBH, QM2;                          \
  396.                                                   \
  397.         __mul_64x64_to_128(ALBH, (A), (B).w[1]);      \
  398.         __mul_64x64_to_128(ALBL, (A), (B).w[0]);       \
  399.                                                   \
  400.         (Ql).w[0] = ALBL.w[0];                          \
  401.     __add_128_64(QM2, ALBH, ALBL.w[1]);           \
  402.         (Ql).w[1] = QM2.w[0];                          \
  403.     Ph = QM2.w[1];                                \
  404. }
  405. #define __mul_64x128_to_192(Q, A, B)              \
  406. {                                                 \
  407. UINT128 ALBL, ALBH, QM2;                          \
  408.                                                   \
  409.         __mul_64x64_to_128(ALBH, (A), (B).w[1]);      \
  410.         __mul_64x64_to_128(ALBL, (A), (B).w[0]);      \
  411.                                                   \
  412.         (Q).w[0] = ALBL.w[0];                         \
  413.     __add_128_64(QM2, ALBH, ALBL.w[1]);           \
  414.         (Q).w[1] = QM2.w[0];                          \
  415.     (Q).w[2] = QM2.w[1];                          \
  416. }
  417. #define __mul_64x128_to192(Q, A, B)          \
  418. {                                             \
  419. UINT128 ALBL, ALBH, QM2;                      \
  420.                                               \
  421.     __mul_64x64_to_128(ALBH, (A), (B).w[1]);  \
  422.     __mul_64x64_to_128(ALBL, (A), (B).w[0]);  \
  423.                                               \
  424.     (Q).w[0] = ALBL.w[0];                    \
  425.     __add_128_64(QM2, ALBH, ALBL.w[1]);       \
  426.     (Q).w[1] = QM2.w[0];                     \
  427.     (Q).w[2] = QM2.w[1];                     \
  428. }
  429. #define __mul_128x128_to_256(P256, A, B)                         \
  430. {                                                                \
  431. UINT128 Qll, Qlh;                                                \
  432. UINT64 Phl, Phh, CY1, CY2;                                         \
  433.                                                                  \
  434.    __mul_64x128_full(Phl, Qll, A.w[0], B);                       \
  435.    __mul_64x128_full(Phh, Qlh, A.w[1], B);                       \
  436.   (P256).w[0] = Qll.w[0];                                        \
  437.            __add_carry_out((P256).w[1],CY1, Qlh.w[0], Qll.w[1]);      \
  438.            __add_carry_in_out((P256).w[2],CY2, Qlh.w[1], Phl, CY1);    \
  439.            (P256).w[3] = Phh + CY2;                                   \
  440. }
  441. //
  442. // For better performance, will check A.w[1] against 0,
  443. //                         but not B.w[1]
  444. // Use this macro accordingly
  445. #define __mul_128x128_to_256_check_A(P256, A, B)                   \
  446. {                                                                  \
  447. UINT128 Qll, Qlh;                                                  \
  448. UINT64 Phl, Phh, CY1, CY2;                                           \
  449.                                                                    \
  450.    __mul_64x128_full(Phl, Qll, A.w[0], B);                          \
  451.   (P256).w[0] = Qll.w[0];                                        \
  452.    if(A.w[1])  {                                                   \
  453.            __mul_64x128_full(Phh, Qlh, A.w[1], B);                     \
  454.            __add_carry_out((P256).w[1],CY1, Qlh.w[0], Qll.w[1]);      \
  455.            __add_carry_in_out((P256).w[2],CY2, Qlh.w[1], Phl, CY1);   \
  456.            (P256).w[3] = Phh + CY2;   }                              \
  457.    else  {                                                         \
  458.            (P256).w[1] = Qll.w[1];                                  \
  459.            (P256).w[2] = Phl;                                       \
  460.            (P256).w[3] = 0;  }                                      \
  461. }
  462. #define __mul_64x192_to_256(lP, lA, lB)                      \
  463. {                                                         \
  464. UINT128 lP0,lP1,lP2;                                      \
  465. UINT64 lC;                                                 \
  466.         __mul_64x64_to_128(lP0, lA, (lB).w[0]);              \
  467.         __mul_64x64_to_128(lP1, lA, (lB).w[1]);              \
  468.         __mul_64x64_to_128(lP2, lA, (lB).w[2]);              \
  469.         (lP).w[0] = lP0.w[0];                                \
  470.         __add_carry_out((lP).w[1],lC,lP1.w[0],lP0.w[1]);      \
  471.         __add_carry_in_out((lP).w[2],lC,lP2.w[0],lP1.w[1],lC); \
  472.         (lP).w[3] = lP2.w[1] + lC;                           \
  473. }
  474. #define __mul_64x256_to_320(P, A, B)                    \
  475. {                                                       \
  476. UINT128 lP0,lP1,lP2,lP3;                                \
  477. UINT64 lC;                                               \
  478.         __mul_64x64_to_128(lP0, A, (B).w[0]);             \
  479.         __mul_64x64_to_128(lP1, A, (B).w[1]);             \
  480.         __mul_64x64_to_128(lP2, A, (B).w[2]);             \
  481.         __mul_64x64_to_128(lP3, A, (B).w[3]);             \
  482.         (P).w[0] = lP0.w[0];                               \
  483.         __add_carry_out((P).w[1],lC,lP1.w[0],lP0.w[1]);      \
  484.         __add_carry_in_out((P).w[2],lC,lP2.w[0],lP1.w[1],lC); \
  485.         __add_carry_in_out((P).w[3],lC,lP3.w[0],lP2.w[1],lC); \
  486.         (P).w[4] = lP3.w[1] + lC;                          \
  487. }
  488. #define __mul_192x192_to_384(P, A, B)                          \
  489. {                                                              \
  490. UINT256 P0,P1,P2;                                              \
  491. UINT64 CY;                                                      \
  492.         __mul_64x192_to_256(P0, (A).w[0], B);                   \
  493.         __mul_64x192_to_256(P1, (A).w[1], B);                   \
  494.         __mul_64x192_to_256(P2, (A).w[2], B);                   \
  495.         (P).w[0] = P0.w[0];                                  \
  496.         __add_carry_out((P).w[1],CY,P1.w[0],P0.w[1]);      \
  497.         __add_carry_in_out((P).w[2],CY,P1.w[1],P0.w[2],CY); \
  498.         __add_carry_in_out((P).w[3],CY,P1.w[2],P0.w[3],CY); \
  499.         (P).w[4] = P1.w[3] + CY;                              \
  500.         __add_carry_out((P).w[2],CY,P2.w[0],(P).w[2]);     \
  501.         __add_carry_in_out((P).w[3],CY,P2.w[1],(P).w[3],CY);   \
  502.         __add_carry_in_out((P).w[4],CY,P2.w[2],(P).w[4],CY);   \
  503.         (P).w[5] = P2.w[3] + CY;                              \
  504. }
  505. #define __mul_64x320_to_384(P, A, B)                    \
  506. {                                                       \
  507. UINT128 lP0,lP1,lP2,lP3,lP4;                            \
  508. UINT64 lC;                                               \
  509.         __mul_64x64_to_128(lP0, A, (B).w[0]);             \
  510.         __mul_64x64_to_128(lP1, A, (B).w[1]);             \
  511.         __mul_64x64_to_128(lP2, A, (B).w[2]);             \
  512.         __mul_64x64_to_128(lP3, A, (B).w[3]);             \
  513.         __mul_64x64_to_128(lP4, A, (B).w[4]);             \
  514.         (P).w[0] = lP0.w[0];                               \
  515.         __add_carry_out((P).w[1],lC,lP1.w[0],lP0.w[1]);      \
  516.         __add_carry_in_out((P).w[2],lC,lP2.w[0],lP1.w[1],lC); \
  517.         __add_carry_in_out((P).w[3],lC,lP3.w[0],lP2.w[1],lC); \
  518.         __add_carry_in_out((P).w[4],lC,lP4.w[0],lP3.w[1],lC); \
  519.         (P).w[5] = lP4.w[1] + lC;                          \
  520. }
  521. // A*A
  522. // Full 128x128-bit product
  523. #define __sqr128_to_256(P256, A)                                 \
  524. {                                                                \
  525. UINT128 Qll, Qlh, Qhh;                                           \
  526. UINT64 TMP_C1, TMP_C2;                                 \
  527.                                                                  \
  528.    __mul_64x64_to_128(Qhh, A.w[1], A.w[1]);                      \
  529.    __mul_64x64_to_128(Qlh, A.w[0], A.w[1]);                      \
  530.    Qhh.w[1] += (Qlh.w[1]>>63);                                   \
  531.    Qlh.w[1] = (Qlh.w[1]+Qlh.w[1])|(Qlh.w[0]>>63);                \
  532.    Qlh.w[0] += Qlh.w[0];                                         \
  533.    __mul_64x64_to_128(Qll, A.w[0], A.w[0]);                      \
  534.                                                                  \
  535.    __add_carry_out((P256).w[1],TMP_C1, Qlh.w[0], Qll.w[1]);      \
  536.    (P256).w[0] = Qll.w[0];                                       \
  537.    __add_carry_in_out((P256).w[2],TMP_C2, Qlh.w[1], Qhh.w[0], TMP_C1);    \
  538.    (P256).w[3] = Qhh.w[1]+TMP_C2;                                         \
  539. }
  540. #define __mul_128x128_to_256_low_high(PQh, PQl, A, B)            \
  541. {                                                                \
  542. UINT128 Qll, Qlh;                                                \
  543. UINT64 Phl, Phh, C1, C2;                                         \
  544.                                                                  \
  545.    __mul_64x128_full(Phl, Qll, A.w[0], B);                       \
  546.    __mul_64x128_full(Phh, Qlh, A.w[1], B);                       \
  547.   (PQl).w[0] = Qll.w[0];                                        \
  548.            __add_carry_out((PQl).w[1],C1, Qlh.w[0], Qll.w[1]);      \
  549.            __add_carry_in_out((PQh).w[0],C2, Qlh.w[1], Phl, C1);    \
  550.            (PQh).w[1] = Phh + C2;                                   \
  551. }
  552. #define __mul_256x256_to_512(P, A, B)                          \
  553. {                                                              \
  554. UINT512 P0,P1,P2,P3;                                           \
  555. UINT64 CY;                                                      \
  556.         __mul_64x256_to_320(P0, (A).w[0], B);                   \
  557.         __mul_64x256_to_320(P1, (A).w[1], B);                   \
  558.         __mul_64x256_to_320(P2, (A).w[2], B);                   \
  559.         __mul_64x256_to_320(P3, (A).w[3], B);                   \
  560.         (P).w[0] = P0.w[0];                                  \
  561.         __add_carry_out((P).w[1],CY,P1.w[0],P0.w[1]);      \
  562.         __add_carry_in_out((P).w[2],CY,P1.w[1],P0.w[2],CY); \
  563.         __add_carry_in_out((P).w[3],CY,P1.w[2],P0.w[3],CY); \
  564.         __add_carry_in_out((P).w[4],CY,P1.w[3],P0.w[4],CY); \
  565.         (P).w[5] = P1.w[4] + CY;                              \
  566.         __add_carry_out((P).w[2],CY,P2.w[0],(P).w[2]);     \
  567.         __add_carry_in_out((P).w[3],CY,P2.w[1],(P).w[3],CY);   \
  568.         __add_carry_in_out((P).w[4],CY,P2.w[2],(P).w[4],CY);   \
  569.         __add_carry_in_out((P).w[5],CY,P2.w[3],(P).w[5],CY);   \
  570.         (P).w[6] = P2.w[4] + CY;                              \
  571.         __add_carry_out((P).w[3],CY,P3.w[0],(P).w[3]);     \
  572.         __add_carry_in_out((P).w[4],CY,P3.w[1],(P).w[4],CY);   \
  573.         __add_carry_in_out((P).w[5],CY,P3.w[2],(P).w[5],CY);   \
  574.         __add_carry_in_out((P).w[6],CY,P3.w[3],(P).w[6],CY);   \
  575.         (P).w[7] = P3.w[4] + CY;                              \
  576. }
  577. #define __mul_192x256_to_448(P, A, B)                          \
  578. {                                                              \
  579. UINT512 P0,P1,P2;                                           \
  580. UINT64 CY;                                                      \
  581.         __mul_64x256_to_320(P0, (A).w[0], B);                   \
  582.         __mul_64x256_to_320(P1, (A).w[1], B);                   \
  583.         __mul_64x256_to_320(P2, (A).w[2], B);                   \
  584.         (P).w[0] = P0.w[0];                                  \
  585.         __add_carry_out((P).w[1],CY,P1.w[0],P0.w[1]);      \
  586.         __add_carry_in_out((P).w[2],CY,P1.w[1],P0.w[2],CY); \
  587.         __add_carry_in_out((P).w[3],CY,P1.w[2],P0.w[3],CY); \
  588.         __add_carry_in_out((P).w[4],CY,P1.w[3],P0.w[4],CY); \
  589.         (P).w[5] = P1.w[4] + CY;                              \
  590.         __add_carry_out((P).w[2],CY,P2.w[0],(P).w[2]);     \
  591.         __add_carry_in_out((P).w[3],CY,P2.w[1],(P).w[3],CY);   \
  592.         __add_carry_in_out((P).w[4],CY,P2.w[2],(P).w[4],CY);   \
  593.         __add_carry_in_out((P).w[5],CY,P2.w[3],(P).w[5],CY);   \
  594.         (P).w[6] = P2.w[4] + CY;                              \
  595. }
  596. #define __mul_320x320_to_640(P, A, B)                          \
  597. {                                                              \
  598. UINT512 P0,P1,P2,P3;                                           \
  599. UINT64 CY;                                                     \
  600.         __mul_256x256_to_512((P), (A), B);                   \
  601.         __mul_64x256_to_320(P1, (A).w[4], B);                   \
  602.         __mul_64x256_to_320(P2, (B).w[4], A);                   \
  603.         __mul_64x64_to_128(P3, (A).w[4], (B).w[4]);               \
  604.         __add_carry_out((P0).w[0],CY,P1.w[0],P2.w[0]);      \
  605.         __add_carry_in_out((P0).w[1],CY,P1.w[1],P2.w[1],CY); \
  606.         __add_carry_in_out((P0).w[2],CY,P1.w[2],P2.w[2],CY); \
  607.         __add_carry_in_out((P0).w[3],CY,P1.w[3],P2.w[3],CY); \
  608.         __add_carry_in_out((P0).w[4],CY,P1.w[4],P2.w[4],CY); \
  609.         P3.w[1] += CY;                                       \
  610.         __add_carry_out((P).w[4],CY,(P).w[4],P0.w[0]);      \
  611.         __add_carry_in_out((P).w[5],CY,(P).w[5],P0.w[1],CY); \
  612.         __add_carry_in_out((P).w[6],CY,(P).w[6],P0.w[2],CY); \
  613.         __add_carry_in_out((P).w[7],CY,(P).w[7],P0.w[3],CY); \
  614.         __add_carry_in_out((P).w[8],CY,P3.w[0],P0.w[4],CY); \
  615.         (P).w[9] = P3.w[1] + CY;                             \
  616. }
  617. #define __mul_384x384_to_768(P, A, B)                          \
  618. {                                                              \
  619. UINT512 P0,P1,P2,P3;                                           \
  620. UINT64 CY;                                                     \
  621.         __mul_320x320_to_640((P), (A), B);                         \
  622.         __mul_64x320_to_384(P1, (A).w[5], B);                   \
  623.         __mul_64x320_to_384(P2, (B).w[5], A);                   \
  624.         __mul_64x64_to_128(P3, (A).w[5], (B).w[5]);               \
  625.         __add_carry_out((P0).w[0],CY,P1.w[0],P2.w[0]);      \
  626.         __add_carry_in_out((P0).w[1],CY,P1.w[1],P2.w[1],CY); \
  627.         __add_carry_in_out((P0).w[2],CY,P1.w[2],P2.w[2],CY); \
  628.         __add_carry_in_out((P0).w[3],CY,P1.w[3],P2.w[3],CY); \
  629.         __add_carry_in_out((P0).w[4],CY,P1.w[4],P2.w[4],CY); \
  630.         __add_carry_in_out((P0).w[5],CY,P1.w[5],P2.w[5],CY); \
  631.         P3.w[1] += CY;                                       \
  632.         __add_carry_out((P).w[5],CY,(P).w[5],P0.w[0]);      \
  633.         __add_carry_in_out((P).w[6],CY,(P).w[6],P0.w[1],CY); \
  634.         __add_carry_in_out((P).w[7],CY,(P).w[7],P0.w[2],CY); \
  635.         __add_carry_in_out((P).w[8],CY,(P).w[8],P0.w[3],CY); \
  636.         __add_carry_in_out((P).w[9],CY,(P).w[9],P0.w[4],CY); \
  637.         __add_carry_in_out((P).w[10],CY,P3.w[0],P0.w[5],CY); \
  638.         (P).w[11] = P3.w[1] + CY;                             \
  639. }
  640. #define __mul_64x128_short(Ql, A, B)              \
  641. {                                                 \
  642. UINT64 ALBH_L;                                    \
  643.                                                   \
  644.         __mul_64x64_to_64(ALBH_L, (A),(B).w[1]);      \
  645.         __mul_64x64_to_128((Ql), (A), (B).w[0]);       \
  646.                                                   \
  647.         (Ql).w[1] += ALBH_L;                          \
  648. }
  649. #define __scale128_10(D,_TMP)                            \
  650. {                                                        \
  651. UINT128 _TMP2,_TMP8;                                     \
  652.           _TMP2.w[1] = (_TMP.w[1]<<1)|(_TMP.w[0]>>63);       \
  653.           _TMP2.w[0] = _TMP.w[0]<<1;                         \
  654.           _TMP8.w[1] = (_TMP.w[1]<<3)|(_TMP.w[0]>>61);       \
  655.           _TMP8.w[0] = _TMP.w[0]<<3;                         \
  656.           __add_128_128(D, _TMP2, _TMP8);                    \
  657. }
  658. // 64x64-bit product
  659. #define __mul_64x64_to_128MACH(P128, CX64, CY64)  \
  660. {                                                  \
  661.   UINT64 CXH,CXL,CYH,CYL,PL,PH,PM,PM2;     \
  662.   CXH = (CX64) >> 32;                              \
  663.   CXL = (UINT32)(CX64);                            \
  664.   CYH = (CY64) >> 32;                              \
  665.   CYL = (UINT32)(CY64);                            \
  666.   PM = CXH*CYL;                                    \
  667.   PH = CXH*CYH;                                    \
  668.   PL = CXL*CYL;                                    \
  669.   PM2 = CXL*CYH;                                   \
  670.   PH += (PM>>32);                                  \
  671.   PM = (UINT64)((UINT32)PM)+PM2+(PL>>32);          \
  672.   (P128).w[1] = PH + (PM>>32);                     \
  673.   (P128).w[0] = (PM<<32)+(UINT32)PL;                \
  674. }
  675. // 64x64-bit product
  676. #define __mul_64x64_to_128HIGH(P64, CX64, CY64)  \
  677. {                                                  \
  678.   UINT64 CXH,CXL,CYH,CYL,PL,PH,PM,PM2;     \
  679.   CXH = (CX64) >> 32;                              \
  680.   CXL = (UINT32)(CX64);                            \
  681.   CYH = (CY64) >> 32;                              \
  682.   CYL = (UINT32)(CY64);                            \
  683.   PM = CXH*CYL;                                    \
  684.   PH = CXH*CYH;                                    \
  685.   PL = CXL*CYL;                                    \
  686.   PM2 = CXL*CYH;                                   \
  687.   PH += (PM>>32);                                  \
  688.   PM = (UINT64)((UINT32)PM)+PM2+(PL>>32);          \
  689.   P64 = PH + (PM>>32);                     \
  690. }
  691. #define __mul_128x64_to_128(Q128, A64, B128)        \
  692. {                                                  \
  693.   UINT64 ALBH_L;                                   \
  694.   ALBH_L = (A64) * (B128).w[1];                    \
  695.   __mul_64x64_to_128MACH((Q128), (A64), (B128).w[0]);   \
  696.   (Q128).w[1] += ALBH_L;                           \
  697. }
  698. // might simplify by calculating just QM2.w[0]
  699. #define __mul_64x128_to_128(Ql, A, B)           \
  700. {                                                 \
  701.   UINT128 ALBL, ALBH, QM2;                        \
  702.   __mul_64x64_to_128(ALBH, (A), (B).w[1]);        \
  703.   __mul_64x64_to_128(ALBL, (A), (B).w[0]);        \
  704.   (Ql).w[0] = ALBL.w[0];                          \
  705.   __add_128_64(QM2, ALBH, ALBL.w[1]);             \
  706.   (Ql).w[1] = QM2.w[0];                           \
  707. }
  708. /*********************************************************************
  709.  *
  710.  *      BID Pack/Unpack Macros
  711.  *
  712.  *********************************************************************/
  713. /////////////////////////////////////////
  714. // BID64 definitions
  715. ////////////////////////////////////////
  716. #define DECIMAL_MAX_EXPON_64  767
  717. #define DECIMAL_EXPONENT_BIAS 398
  718. #define MAX_FORMAT_DIGITS     16
  719. /////////////////////////////////////////
  720. // BID128 definitions
  721. ////////////////////////////////////////
  722. #define DECIMAL_MAX_EXPON_128  12287
  723. #define DECIMAL_EXPONENT_BIAS_128  6176
  724. #define MAX_FORMAT_DIGITS_128      34
  725. /////////////////////////////////////////
  726. // BID32 definitions
  727. ////////////////////////////////////////
  728. #define DECIMAL_MAX_EXPON_32  191
  729. #define DECIMAL_EXPONENT_BIAS_32  101
  730. #define MAX_FORMAT_DIGITS_32      7
  731. ////////////////////////////////////////
  732. // Constant Definitions
  733. ///////////////////////////////////////
  734. #define SPECIAL_ENCODING_MASK64 0x6000000000000000ull
  735. #define INFINITY_MASK64         0x7800000000000000ull
  736. #define SINFINITY_MASK64        0xf800000000000000ull
  737. #define SSNAN_MASK64            0xfc00000000000000ull
  738. #define NAN_MASK64              0x7c00000000000000ull
  739. #define SNAN_MASK64             0x7e00000000000000ull
  740. #define QUIET_MASK64            0xfdffffffffffffffull
  741. #define LARGE_COEFF_MASK64      0x0007ffffffffffffull
  742. #define LARGE_COEFF_HIGH_BIT64  0x0020000000000000ull
  743. #define SMALL_COEFF_MASK64      0x001fffffffffffffull
  744. #define EXPONENT_MASK64         0x3ff
  745. #define EXPONENT_SHIFT_LARGE64  51
  746. #define EXPONENT_SHIFT_SMALL64  53
  747. #define LARGEST_BID64           0x77fb86f26fc0ffffull
  748. #define SMALLEST_BID64          0xf7fb86f26fc0ffffull
  749. #define SMALL_COEFF_MASK128     0x0001ffffffffffffull
  750. #define LARGE_COEFF_MASK128     0x00007fffffffffffull
  751. #define EXPONENT_MASK128        0x3fff
  752. #define LARGEST_BID128_HIGH     0x5fffed09bead87c0ull
  753. #define LARGEST_BID128_LOW      0x378d8e63ffffffffull
  754. #define SPECIAL_ENCODING_MASK32 0x60000000ul
  755. #define INFINITY_MASK32         0x78000000ul
  756. #define LARGE_COEFF_MASK32      0x007ffffful
  757. #define LARGE_COEFF_HIGH_BIT32  0x00800000ul
  758. #define SMALL_COEFF_MASK32      0x001ffffful
  759. #define EXPONENT_MASK32         0xff
  760. #define LARGEST_BID32           0x77f8967f
  761. #define NAN_MASK32              0x7c000000
  762. #define SNAN_MASK32             0x7e000000
  763. #define MASK_BINARY_EXPONENT  0x7ff0000000000000ull
  764. #define BINARY_EXPONENT_BIAS  0x3ff
  765. #define UPPER_EXPON_LIMIT     51
  766. // data needed for BID pack/unpack macros
  767. extern UINT64 round_const_table[][19];
  768. extern UINT128 reciprocals10_128[];
  769. extern int recip_scale[];
  770. extern UINT128 power10_table_128[];
  771. extern int estimate_decimal_digits[];
  772. extern int estimate_bin_expon[];
  773. extern UINT64 power10_index_binexp[];
  774. extern int short_recip_scale[];
  775. extern UINT64 reciprocals10_64[];
  776. extern UINT128 power10_index_binexp_128[];
  777. extern UINT128 round_const_table_128[][36];
  778.  
  779.  
  780. //////////////////////////////////////////////
  781. //  Status Flag Handling
  782. /////////////////////////////////////////////
  783. #define __set_status_flags(fpsc, status)  *(fpsc) |= status
  784. #define is_inexact(fpsc)  ((*(fpsc))&INEXACT_EXCEPTION)
  785.  
  786. __BID_INLINE__ UINT64
  787. unpack_BID64 (UINT64 * psign_x, int *pexponent_x,
  788.               UINT64 * pcoefficient_x, UINT64 x) {
  789.   UINT64 tmp, coeff;
  790.  
  791.   *psign_x = x & 0x8000000000000000ull;
  792.  
  793.   if ((x & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64) {
  794.     // special encodings
  795.     // coefficient
  796.     coeff = (x & LARGE_COEFF_MASK64) | LARGE_COEFF_HIGH_BIT64;
  797.  
  798.     if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
  799.       *pexponent_x = 0;
  800.       *pcoefficient_x = x & 0xfe03ffffffffffffull;
  801.       if ((x & 0x0003ffffffffffffull) >= 1000000000000000ull)
  802.         *pcoefficient_x = x & 0xfe00000000000000ull;
  803.       if ((x & NAN_MASK64) == INFINITY_MASK64)
  804.         *pcoefficient_x = x & SINFINITY_MASK64;
  805.       return 0; // NaN or Infinity
  806.     }
  807.     // check for non-canonical values
  808.     if (coeff >= 10000000000000000ull)
  809.       coeff = 0;
  810.     *pcoefficient_x = coeff;
  811.     // get exponent
  812.     tmp = x >> EXPONENT_SHIFT_LARGE64;
  813.     *pexponent_x = (int) (tmp & EXPONENT_MASK64);
  814.     return coeff;
  815.   }
  816.   // exponent
  817.   tmp = x >> EXPONENT_SHIFT_SMALL64;
  818.   *pexponent_x = (int) (tmp & EXPONENT_MASK64);
  819.   // coefficient
  820.   *pcoefficient_x = (x & SMALL_COEFF_MASK64);
  821.  
  822.   return *pcoefficient_x;
  823. }
  824.  
  825. //
  826. //   BID64 pack macro (general form)
  827. //
  828. __BID_INLINE__ UINT64
  829. get_BID64 (UINT64 sgn, int expon, UINT64 coeff, int rmode,
  830.            unsigned *fpsc) {
  831.   UINT128 Stemp, Q_low;
  832.   UINT64 QH, r, mask, C64, remainder_h, CY, carry;
  833.   int extra_digits, amount, amount2;
  834.   unsigned status;
  835.  
  836.   if (coeff > 9999999999999999ull) {
  837.     expon++;
  838.     coeff = 1000000000000000ull;
  839.   }
  840.   // check for possible underflow/overflow
  841.   if (((unsigned) expon) >= 3 * 256) {
  842.     if (expon < 0) {
  843.       // underflow
  844.       if (expon + MAX_FORMAT_DIGITS < 0) {
  845. #ifdef SET_STATUS_FLAGS
  846.         __set_status_flags (fpsc,
  847.                             UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  848. #endif
  849. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  850. #ifndef IEEE_ROUND_NEAREST
  851.         if (rmode == ROUNDING_DOWN && sgn)
  852.           return 0x8000000000000001ull;
  853.         if (rmode == ROUNDING_UP && !sgn)
  854.           return 1ull;
  855. #endif
  856. #endif
  857.         // result is 0
  858.         return sgn;
  859.       }
  860. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  861. #ifndef IEEE_ROUND_NEAREST
  862.       if (sgn && (unsigned) (rmode - 1) < 2)
  863.         rmode = 3 - rmode;
  864. #endif
  865. #endif
  866.       // get digits to be shifted out
  867.       extra_digits = -expon;
  868.       coeff += round_const_table[rmode][extra_digits];
  869.  
  870.       // get coeff*(2^M[extra_digits])/10^extra_digits
  871.       __mul_64x128_full (QH, Q_low, coeff,
  872.                          reciprocals10_128[extra_digits]);
  873.  
  874.       // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
  875.       amount = recip_scale[extra_digits];
  876.  
  877.       C64 = QH >> amount;
  878.  
  879. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  880. #ifndef IEEE_ROUND_NEAREST
  881.       if (rmode == 0)   //ROUNDING_TO_NEAREST
  882. #endif
  883.         if (C64 & 1) {
  884.           // check whether fractional part of initial_P/10^extra_digits is exactly .5
  885.  
  886.           // get remainder
  887.           amount2 = 64 - amount;
  888.           remainder_h = 0;
  889.           remainder_h--;
  890.           remainder_h >>= amount2;
  891.           remainder_h = remainder_h & QH;
  892.  
  893.           if (!remainder_h
  894.               && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
  895.                   || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
  896.                       && Q_low.w[0] <
  897.                       reciprocals10_128[extra_digits].w[0]))) {
  898.             C64--;
  899.           }
  900.         }
  901. #endif
  902.  
  903. #ifdef SET_STATUS_FLAGS
  904.  
  905.       if (is_inexact (fpsc))
  906.         __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
  907.       else {
  908.         status = INEXACT_EXCEPTION;
  909.         // get remainder
  910.         remainder_h = QH << (64 - amount);
  911.  
  912.         switch (rmode) {
  913.         case ROUNDING_TO_NEAREST:
  914.         case ROUNDING_TIES_AWAY:
  915.           // test whether fractional part is 0
  916.           if (remainder_h == 0x8000000000000000ull
  917.               && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
  918.                   || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
  919.                       && Q_low.w[0] <
  920.                       reciprocals10_128[extra_digits].w[0])))
  921.             status = EXACT_STATUS;
  922.           break;
  923.         case ROUNDING_DOWN:
  924.         case ROUNDING_TO_ZERO:
  925.           if (!remainder_h
  926.               && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
  927.                   || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
  928.                       && Q_low.w[0] <
  929.                       reciprocals10_128[extra_digits].w[0])))
  930.             status = EXACT_STATUS;
  931.           break;
  932.         default:
  933.           // round up
  934.           __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
  935.                            reciprocals10_128[extra_digits].w[0]);
  936.           __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
  937.                               reciprocals10_128[extra_digits].w[1], CY);
  938.           if ((remainder_h >> (64 - amount)) + carry >=
  939.               (((UINT64) 1) << amount))
  940.             status = EXACT_STATUS;
  941.         }
  942.  
  943.         if (status != EXACT_STATUS)
  944.           __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status);
  945.       }
  946.  
  947. #endif
  948.  
  949.       return sgn | C64;
  950.     }
  951.     while (coeff < 1000000000000000ull && expon >= 3 * 256) {
  952.       expon--;
  953.       coeff = (coeff << 3) + (coeff << 1);
  954.     }
  955.     if (expon > DECIMAL_MAX_EXPON_64) {
  956. #ifdef SET_STATUS_FLAGS
  957.       __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  958. #endif
  959.       // overflow
  960.       r = sgn | INFINITY_MASK64;
  961.       switch (rmode) {
  962.       case ROUNDING_DOWN:
  963.         if (!sgn)
  964.           r = LARGEST_BID64;
  965.         break;
  966.       case ROUNDING_TO_ZERO:
  967.         r = sgn | LARGEST_BID64;
  968.         break;
  969.       case ROUNDING_UP:
  970.         // round up
  971.         if (sgn)
  972.           r = SMALLEST_BID64;
  973.       }
  974.       return r;
  975.     }
  976.   }
  977.  
  978.   mask = 1;
  979.   mask <<= EXPONENT_SHIFT_SMALL64;
  980.  
  981.   // check whether coefficient fits in 10*5+3 bits
  982.   if (coeff < mask) {
  983.     r = expon;
  984.     r <<= EXPONENT_SHIFT_SMALL64;
  985.     r |= (coeff | sgn);
  986.     return r;
  987.   }
  988.   // special format
  989.  
  990.   // eliminate the case coeff==10^16 after rounding
  991.   if (coeff == 10000000000000000ull) {
  992.     r = expon + 1;
  993.     r <<= EXPONENT_SHIFT_SMALL64;
  994.     r |= (1000000000000000ull | sgn);
  995.     return r;
  996.   }
  997.  
  998.   r = expon;
  999.   r <<= EXPONENT_SHIFT_LARGE64;
  1000.   r |= (sgn | SPECIAL_ENCODING_MASK64);
  1001.   // add coeff, without leading bits
  1002.   mask = (mask >> 2) - 1;
  1003.   coeff &= mask;
  1004.   r |= coeff;
  1005.  
  1006.   return r;
  1007. }
  1008.  
  1009.  
  1010.  
  1011.  
  1012. //
  1013. //   No overflow/underflow checking
  1014. //
  1015. __BID_INLINE__ UINT64
  1016. fast_get_BID64 (UINT64 sgn, int expon, UINT64 coeff) {
  1017.   UINT64 r, mask;
  1018.  
  1019.   mask = 1;
  1020.   mask <<= EXPONENT_SHIFT_SMALL64;
  1021.  
  1022.   // check whether coefficient fits in 10*5+3 bits
  1023.   if (coeff < mask) {
  1024.     r = expon;
  1025.     r <<= EXPONENT_SHIFT_SMALL64;
  1026.     r |= (coeff | sgn);
  1027.     return r;
  1028.   }
  1029.   // special format
  1030.  
  1031.   // eliminate the case coeff==10^16 after rounding
  1032.   if (coeff == 10000000000000000ull) {
  1033.     r = expon + 1;
  1034.     r <<= EXPONENT_SHIFT_SMALL64;
  1035.     r |= (1000000000000000ull | sgn);
  1036.     return r;
  1037.   }
  1038.  
  1039.   r = expon;
  1040.   r <<= EXPONENT_SHIFT_LARGE64;
  1041.   r |= (sgn | SPECIAL_ENCODING_MASK64);
  1042.   // add coeff, without leading bits
  1043.   mask = (mask >> 2) - 1;
  1044.   coeff &= mask;
  1045.   r |= coeff;
  1046.  
  1047.   return r;
  1048. }
  1049.  
  1050.  
  1051. //
  1052. //   no underflow checking
  1053. //
  1054. __BID_INLINE__ UINT64
  1055. fast_get_BID64_check_OF (UINT64 sgn, int expon, UINT64 coeff, int rmode,
  1056.                          unsigned *fpsc) {
  1057.   UINT64 r, mask;
  1058.  
  1059.   if (((unsigned) expon) >= 3 * 256 - 1) {
  1060.     if ((expon == 3 * 256 - 1) && coeff == 10000000000000000ull) {
  1061.       expon = 3 * 256;
  1062.       coeff = 1000000000000000ull;
  1063.     }
  1064.  
  1065.     if (((unsigned) expon) >= 3 * 256) {
  1066.       while (coeff < 1000000000000000ull && expon >= 3 * 256) {
  1067.         expon--;
  1068.         coeff = (coeff << 3) + (coeff << 1);
  1069.       }
  1070.       if (expon > DECIMAL_MAX_EXPON_64) {
  1071. #ifdef SET_STATUS_FLAGS
  1072.         __set_status_flags (fpsc,
  1073.                             OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  1074. #endif
  1075.         // overflow
  1076.         r = sgn | INFINITY_MASK64;
  1077.         switch (rmode) {
  1078.         case ROUNDING_DOWN:
  1079.           if (!sgn)
  1080.             r = LARGEST_BID64;
  1081.           break;
  1082.         case ROUNDING_TO_ZERO:
  1083.           r = sgn | LARGEST_BID64;
  1084.           break;
  1085.         case ROUNDING_UP:
  1086.           // round up
  1087.           if (sgn)
  1088.             r = SMALLEST_BID64;
  1089.         }
  1090.         return r;
  1091.       }
  1092.     }
  1093.   }
  1094.  
  1095.   mask = 1;
  1096.   mask <<= EXPONENT_SHIFT_SMALL64;
  1097.  
  1098.   // check whether coefficient fits in 10*5+3 bits
  1099.   if (coeff < mask) {
  1100.     r = expon;
  1101.     r <<= EXPONENT_SHIFT_SMALL64;
  1102.     r |= (coeff | sgn);
  1103.     return r;
  1104.   }
  1105.   // special format
  1106.  
  1107.   // eliminate the case coeff==10^16 after rounding
  1108.   if (coeff == 10000000000000000ull) {
  1109.     r = expon + 1;
  1110.     r <<= EXPONENT_SHIFT_SMALL64;
  1111.     r |= (1000000000000000ull | sgn);
  1112.     return r;
  1113.   }
  1114.  
  1115.   r = expon;
  1116.   r <<= EXPONENT_SHIFT_LARGE64;
  1117.   r |= (sgn | SPECIAL_ENCODING_MASK64);
  1118.   // add coeff, without leading bits
  1119.   mask = (mask >> 2) - 1;
  1120.   coeff &= mask;
  1121.   r |= coeff;
  1122.  
  1123.   return r;
  1124. }
  1125.  
  1126.  
  1127. //
  1128. //   No overflow/underflow checking
  1129. //   or checking for coefficients equal to 10^16 (after rounding)
  1130. //
  1131. __BID_INLINE__ UINT64
  1132. very_fast_get_BID64 (UINT64 sgn, int expon, UINT64 coeff) {
  1133.   UINT64 r, mask;
  1134.  
  1135.   mask = 1;
  1136.   mask <<= EXPONENT_SHIFT_SMALL64;
  1137.  
  1138.   // check whether coefficient fits in 10*5+3 bits
  1139.   if (coeff < mask) {
  1140.     r = expon;
  1141.     r <<= EXPONENT_SHIFT_SMALL64;
  1142.     r |= (coeff | sgn);
  1143.     return r;
  1144.   }
  1145.   // special format
  1146.   r = expon;
  1147.   r <<= EXPONENT_SHIFT_LARGE64;
  1148.   r |= (sgn | SPECIAL_ENCODING_MASK64);
  1149.   // add coeff, without leading bits
  1150.   mask = (mask >> 2) - 1;
  1151.   coeff &= mask;
  1152.   r |= coeff;
  1153.  
  1154.   return r;
  1155. }
  1156.  
  1157. //
  1158. //   No overflow/underflow checking or checking for coefficients above 2^53
  1159. //
  1160. __BID_INLINE__ UINT64
  1161. very_fast_get_BID64_small_mantissa (UINT64 sgn, int expon, UINT64 coeff) {
  1162.   // no UF/OF
  1163.   UINT64 r;
  1164.  
  1165.   r = expon;
  1166.   r <<= EXPONENT_SHIFT_SMALL64;
  1167.   r |= (coeff | sgn);
  1168.   return r;
  1169. }
  1170.  
  1171.  
  1172. //
  1173. // This pack macro is used when underflow is known to occur
  1174. //
  1175. __BID_INLINE__ UINT64
  1176. get_BID64_UF (UINT64 sgn, int expon, UINT64 coeff, UINT64 R, int rmode,
  1177.               unsigned *fpsc) {
  1178.   UINT128 C128, Q_low, Stemp;
  1179.   UINT64 C64, remainder_h, QH, carry, CY;
  1180.   int extra_digits, amount, amount2;
  1181.   unsigned status;
  1182.  
  1183.   // underflow
  1184.   if (expon + MAX_FORMAT_DIGITS < 0) {
  1185. #ifdef SET_STATUS_FLAGS
  1186.     __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  1187. #endif
  1188. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1189. #ifndef IEEE_ROUND_NEAREST
  1190.     if (rmode == ROUNDING_DOWN && sgn)
  1191.       return 0x8000000000000001ull;
  1192.     if (rmode == ROUNDING_UP && !sgn)
  1193.       return 1ull;
  1194. #endif
  1195. #endif
  1196.     // result is 0
  1197.     return sgn;
  1198.   }
  1199.   // 10*coeff
  1200.   coeff = (coeff << 3) + (coeff << 1);
  1201. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1202. #ifndef IEEE_ROUND_NEAREST
  1203.   if (sgn && (unsigned) (rmode - 1) < 2)
  1204.     rmode = 3 - rmode;
  1205. #endif
  1206. #endif
  1207.   if (R)
  1208.     coeff |= 1;
  1209.   // get digits to be shifted out
  1210.   extra_digits = 1 - expon;
  1211.   C128.w[0] = coeff + round_const_table[rmode][extra_digits];
  1212.  
  1213.   // get coeff*(2^M[extra_digits])/10^extra_digits
  1214.   __mul_64x128_full (QH, Q_low, C128.w[0],
  1215.                      reciprocals10_128[extra_digits]);
  1216.  
  1217.   // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
  1218.   amount = recip_scale[extra_digits];
  1219.  
  1220.   C64 = QH >> amount;
  1221.   //__shr_128(C128, Q_high, amount);
  1222.  
  1223. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1224. #ifndef IEEE_ROUND_NEAREST
  1225.   if (rmode == 0)       //ROUNDING_TO_NEAREST
  1226. #endif
  1227.     if (C64 & 1) {
  1228.       // check whether fractional part of initial_P/10^extra_digits is exactly .5
  1229.  
  1230.       // get remainder
  1231.       amount2 = 64 - amount;
  1232.       remainder_h = 0;
  1233.       remainder_h--;
  1234.       remainder_h >>= amount2;
  1235.       remainder_h = remainder_h & QH;
  1236.  
  1237.       if (!remainder_h
  1238.           && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
  1239.               || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
  1240.                   && Q_low.w[0] <
  1241.                   reciprocals10_128[extra_digits].w[0]))) {
  1242.         C64--;
  1243.       }
  1244.     }
  1245. #endif
  1246.  
  1247. #ifdef SET_STATUS_FLAGS
  1248.  
  1249.   if (is_inexact (fpsc))
  1250.     __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
  1251.   else {
  1252.     status = INEXACT_EXCEPTION;
  1253.     // get remainder
  1254.     remainder_h = QH << (64 - amount);
  1255.  
  1256.     switch (rmode) {
  1257.     case ROUNDING_TO_NEAREST:
  1258.     case ROUNDING_TIES_AWAY:
  1259.       // test whether fractional part is 0
  1260.       if (remainder_h == 0x8000000000000000ull
  1261.           && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
  1262.               || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
  1263.                   && Q_low.w[0] <
  1264.                   reciprocals10_128[extra_digits].w[0])))
  1265.         status = EXACT_STATUS;
  1266.       break;
  1267.     case ROUNDING_DOWN:
  1268.     case ROUNDING_TO_ZERO:
  1269.       if (!remainder_h
  1270.           && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
  1271.               || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
  1272.                   && Q_low.w[0] <
  1273.                   reciprocals10_128[extra_digits].w[0])))
  1274.         status = EXACT_STATUS;
  1275.       break;
  1276.     default:
  1277.       // round up
  1278.       __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
  1279.                        reciprocals10_128[extra_digits].w[0]);
  1280.       __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
  1281.                           reciprocals10_128[extra_digits].w[1], CY);
  1282.       if ((remainder_h >> (64 - amount)) + carry >=
  1283.           (((UINT64) 1) << amount))
  1284.         status = EXACT_STATUS;
  1285.     }
  1286.  
  1287.     if (status != EXACT_STATUS)
  1288.       __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status);
  1289.   }
  1290.  
  1291. #endif
  1292.  
  1293.   return sgn | C64;
  1294.  
  1295. }
  1296.  
  1297.  
  1298.  
  1299. //
  1300. //   This pack macro doesnot check for coefficients above 2^53
  1301. //
  1302. __BID_INLINE__ UINT64
  1303. get_BID64_small_mantissa (UINT64 sgn, int expon, UINT64 coeff,
  1304.                           int rmode, unsigned *fpsc) {
  1305.   UINT128 C128, Q_low, Stemp;
  1306.   UINT64 r, mask, C64, remainder_h, QH, carry, CY;
  1307.   int extra_digits, amount, amount2;
  1308.   unsigned status;
  1309.  
  1310.   // check for possible underflow/overflow
  1311.   if (((unsigned) expon) >= 3 * 256) {
  1312.     if (expon < 0) {
  1313.       // underflow
  1314.       if (expon + MAX_FORMAT_DIGITS < 0) {
  1315. #ifdef SET_STATUS_FLAGS
  1316.         __set_status_flags (fpsc,
  1317.                             UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  1318. #endif
  1319. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1320. #ifndef IEEE_ROUND_NEAREST
  1321.         if (rmode == ROUNDING_DOWN && sgn)
  1322.           return 0x8000000000000001ull;
  1323.         if (rmode == ROUNDING_UP && !sgn)
  1324.           return 1ull;
  1325. #endif
  1326. #endif
  1327.         // result is 0
  1328.         return sgn;
  1329.       }
  1330. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1331. #ifndef IEEE_ROUND_NEAREST
  1332.       if (sgn && (unsigned) (rmode - 1) < 2)
  1333.         rmode = 3 - rmode;
  1334. #endif
  1335. #endif
  1336.       // get digits to be shifted out
  1337.       extra_digits = -expon;
  1338.       C128.w[0] = coeff + round_const_table[rmode][extra_digits];
  1339.  
  1340.       // get coeff*(2^M[extra_digits])/10^extra_digits
  1341.       __mul_64x128_full (QH, Q_low, C128.w[0],
  1342.                          reciprocals10_128[extra_digits]);
  1343.  
  1344.       // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
  1345.       amount = recip_scale[extra_digits];
  1346.  
  1347.       C64 = QH >> amount;
  1348.  
  1349. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1350. #ifndef IEEE_ROUND_NEAREST
  1351.       if (rmode == 0)   //ROUNDING_TO_NEAREST
  1352. #endif
  1353.         if (C64 & 1) {
  1354.           // check whether fractional part of initial_P/10^extra_digits is exactly .5
  1355.  
  1356.           // get remainder
  1357.           amount2 = 64 - amount;
  1358.           remainder_h = 0;
  1359.           remainder_h--;
  1360.           remainder_h >>= amount2;
  1361.           remainder_h = remainder_h & QH;
  1362.  
  1363.           if (!remainder_h
  1364.               && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
  1365.                   || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
  1366.                       && Q_low.w[0] <
  1367.                       reciprocals10_128[extra_digits].w[0]))) {
  1368.             C64--;
  1369.           }
  1370.         }
  1371. #endif
  1372.  
  1373. #ifdef SET_STATUS_FLAGS
  1374.  
  1375.       if (is_inexact (fpsc))
  1376.         __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
  1377.       else {
  1378.         status = INEXACT_EXCEPTION;
  1379.         // get remainder
  1380.         remainder_h = QH << (64 - amount);
  1381.  
  1382.         switch (rmode) {
  1383.         case ROUNDING_TO_NEAREST:
  1384.         case ROUNDING_TIES_AWAY:
  1385.           // test whether fractional part is 0
  1386.           if (remainder_h == 0x8000000000000000ull
  1387.               && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
  1388.                   || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
  1389.                       && Q_low.w[0] <
  1390.                       reciprocals10_128[extra_digits].w[0])))
  1391.             status = EXACT_STATUS;
  1392.           break;
  1393.         case ROUNDING_DOWN:
  1394.         case ROUNDING_TO_ZERO:
  1395.           if (!remainder_h
  1396.               && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
  1397.                   || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
  1398.                       && Q_low.w[0] <
  1399.                       reciprocals10_128[extra_digits].w[0])))
  1400.             status = EXACT_STATUS;
  1401.           break;
  1402.         default:
  1403.           // round up
  1404.           __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
  1405.                            reciprocals10_128[extra_digits].w[0]);
  1406.           __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
  1407.                               reciprocals10_128[extra_digits].w[1], CY);
  1408.           if ((remainder_h >> (64 - amount)) + carry >=
  1409.               (((UINT64) 1) << amount))
  1410.             status = EXACT_STATUS;
  1411.         }
  1412.  
  1413.         if (status != EXACT_STATUS)
  1414.           __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status);
  1415.       }
  1416.  
  1417. #endif
  1418.  
  1419.       return sgn | C64;
  1420.     }
  1421.  
  1422.     while (coeff < 1000000000000000ull && expon >= 3 * 256) {
  1423.       expon--;
  1424.       coeff = (coeff << 3) + (coeff << 1);
  1425.     }
  1426.     if (expon > DECIMAL_MAX_EXPON_64) {
  1427. #ifdef SET_STATUS_FLAGS
  1428.       __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  1429. #endif
  1430.       // overflow
  1431.       r = sgn | INFINITY_MASK64;
  1432.       switch (rmode) {
  1433.       case ROUNDING_DOWN:
  1434.         if (!sgn)
  1435.           r = LARGEST_BID64;
  1436.         break;
  1437.       case ROUNDING_TO_ZERO:
  1438.         r = sgn | LARGEST_BID64;
  1439.         break;
  1440.       case ROUNDING_UP:
  1441.         // round up
  1442.         if (sgn)
  1443.           r = SMALLEST_BID64;
  1444.       }
  1445.       return r;
  1446.     } else {
  1447.       mask = 1;
  1448.       mask <<= EXPONENT_SHIFT_SMALL64;
  1449.       if (coeff >= mask) {
  1450.         r = expon;
  1451.         r <<= EXPONENT_SHIFT_LARGE64;
  1452.         r |= (sgn | SPECIAL_ENCODING_MASK64);
  1453.         // add coeff, without leading bits
  1454.         mask = (mask >> 2) - 1;
  1455.         coeff &= mask;
  1456.         r |= coeff;
  1457.         return r;
  1458.       }
  1459.     }
  1460.   }
  1461.  
  1462.   r = expon;
  1463.   r <<= EXPONENT_SHIFT_SMALL64;
  1464.   r |= (coeff | sgn);
  1465.  
  1466.   return r;
  1467. }
  1468.  
  1469.  
  1470. /*****************************************************************************
  1471. *
  1472. *    BID128 pack/unpack macros
  1473. *
  1474. *****************************************************************************/
  1475.  
  1476. //
  1477. //   Macro for handling BID128 underflow
  1478. //         sticky bit given as additional argument
  1479. //
  1480. __BID_INLINE__ UINT128 *
  1481. handle_UF_128_rem (UINT128 * pres, UINT64 sgn, int expon, UINT128 CQ,
  1482.                    UINT64 R, unsigned *prounding_mode, unsigned *fpsc) {
  1483.   UINT128 T128, TP128, Qh, Ql, Qh1, Stemp, Tmp, Tmp1, CQ2, CQ8;
  1484.   UINT64 carry, CY;
  1485.   int ed2, amount;
  1486.   unsigned rmode, status;
  1487.  
  1488.   // UF occurs
  1489.   if (expon + MAX_FORMAT_DIGITS_128 < 0) {
  1490. #ifdef SET_STATUS_FLAGS
  1491.     __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  1492. #endif
  1493.     pres->w[1] = sgn;
  1494.     pres->w[0] = 0;
  1495. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1496. #ifndef IEEE_ROUND_NEAREST
  1497.     if ((sgn && *prounding_mode == ROUNDING_DOWN)
  1498.         || (!sgn && *prounding_mode == ROUNDING_UP))
  1499.       pres->w[0] = 1ull;
  1500. #endif
  1501. #endif
  1502.     return pres;
  1503.   }
  1504.   // CQ *= 10
  1505.   CQ2.w[1] = (CQ.w[1] << 1) | (CQ.w[0] >> 63);
  1506.   CQ2.w[0] = CQ.w[0] << 1;
  1507.   CQ8.w[1] = (CQ.w[1] << 3) | (CQ.w[0] >> 61);
  1508.   CQ8.w[0] = CQ.w[0] << 3;
  1509.   __add_128_128 (CQ, CQ2, CQ8);
  1510.  
  1511.   // add remainder
  1512.   if (R)
  1513.     CQ.w[0] |= 1;
  1514.  
  1515.   ed2 = 1 - expon;
  1516.   // add rounding constant to CQ
  1517. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1518. #ifndef IEEE_ROUND_NEAREST
  1519.   rmode = *prounding_mode;
  1520.   if (sgn && (unsigned) (rmode - 1) < 2)
  1521.     rmode = 3 - rmode;
  1522. #else
  1523.   rmode = 0;
  1524. #endif
  1525. #else
  1526.   rmode = 0;
  1527. #endif
  1528.   T128 = round_const_table_128[rmode][ed2];
  1529.   __add_carry_out (CQ.w[0], carry, T128.w[0], CQ.w[0]);
  1530.   CQ.w[1] = CQ.w[1] + T128.w[1] + carry;
  1531.  
  1532.   TP128 = reciprocals10_128[ed2];
  1533.   __mul_128x128_full (Qh, Ql, CQ, TP128);
  1534.   amount = recip_scale[ed2];
  1535.  
  1536.   if (amount >= 64) {
  1537.     CQ.w[0] = Qh.w[1] >> (amount - 64);
  1538.     CQ.w[1] = 0;
  1539.   } else {
  1540.     __shr_128 (CQ, Qh, amount);
  1541.   }
  1542.  
  1543.   expon = 0;
  1544.  
  1545. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1546. #ifndef IEEE_ROUND_NEAREST
  1547.   if (!(*prounding_mode))
  1548. #endif
  1549.     if (CQ.w[0] & 1) {
  1550.       // check whether fractional part of initial_P/10^ed1 is exactly .5
  1551.  
  1552.       // get remainder
  1553.       __shl_128_long (Qh1, Qh, (128 - amount));
  1554.  
  1555.       if (!Qh1.w[1] && !Qh1.w[0]
  1556.           && (Ql.w[1] < reciprocals10_128[ed2].w[1]
  1557.               || (Ql.w[1] == reciprocals10_128[ed2].w[1]
  1558.                   && Ql.w[0] < reciprocals10_128[ed2].w[0]))) {
  1559.         CQ.w[0]--;
  1560.       }
  1561.     }
  1562. #endif
  1563.  
  1564. #ifdef SET_STATUS_FLAGS
  1565.  
  1566.   if (is_inexact (fpsc))
  1567.     __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
  1568.   else {
  1569.     status = INEXACT_EXCEPTION;
  1570.     // get remainder
  1571.     __shl_128_long (Qh1, Qh, (128 - amount));
  1572.  
  1573.     switch (rmode) {
  1574.     case ROUNDING_TO_NEAREST:
  1575.     case ROUNDING_TIES_AWAY:
  1576.       // test whether fractional part is 0
  1577.       if (Qh1.w[1] == 0x8000000000000000ull && (!Qh1.w[0])
  1578.           && (Ql.w[1] < reciprocals10_128[ed2].w[1]
  1579.               || (Ql.w[1] == reciprocals10_128[ed2].w[1]
  1580.                   && Ql.w[0] < reciprocals10_128[ed2].w[0])))
  1581.         status = EXACT_STATUS;
  1582.       break;
  1583.     case ROUNDING_DOWN:
  1584.     case ROUNDING_TO_ZERO:
  1585.       if ((!Qh1.w[1]) && (!Qh1.w[0])
  1586.           && (Ql.w[1] < reciprocals10_128[ed2].w[1]
  1587.               || (Ql.w[1] == reciprocals10_128[ed2].w[1]
  1588.                   && Ql.w[0] < reciprocals10_128[ed2].w[0])))
  1589.         status = EXACT_STATUS;
  1590.       break;
  1591.     default:
  1592.       // round up
  1593.       __add_carry_out (Stemp.w[0], CY, Ql.w[0],
  1594.                        reciprocals10_128[ed2].w[0]);
  1595.       __add_carry_in_out (Stemp.w[1], carry, Ql.w[1],
  1596.                           reciprocals10_128[ed2].w[1], CY);
  1597.       __shr_128_long (Qh, Qh1, (128 - amount));
  1598.       Tmp.w[0] = 1;
  1599.       Tmp.w[1] = 0;
  1600.       __shl_128_long (Tmp1, Tmp, amount);
  1601.       Qh.w[0] += carry;
  1602.       if (Qh.w[0] < carry)
  1603.         Qh.w[1]++;
  1604.       if (__unsigned_compare_ge_128 (Qh, Tmp1))
  1605.         status = EXACT_STATUS;
  1606.     }
  1607.  
  1608.     if (status != EXACT_STATUS)
  1609.       __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status);
  1610.   }
  1611.  
  1612. #endif
  1613.  
  1614.   pres->w[1] = sgn | CQ.w[1];
  1615.   pres->w[0] = CQ.w[0];
  1616.  
  1617.   return pres;
  1618.  
  1619. }
  1620.  
  1621.  
  1622. //
  1623. //   Macro for handling BID128 underflow
  1624. //
  1625. __BID_INLINE__ UINT128 *
  1626. handle_UF_128 (UINT128 * pres, UINT64 sgn, int expon, UINT128 CQ,
  1627.                unsigned *prounding_mode, unsigned *fpsc) {
  1628.   UINT128 T128, TP128, Qh, Ql, Qh1, Stemp, Tmp, Tmp1;
  1629.   UINT64 carry, CY;
  1630.   int ed2, amount;
  1631.   unsigned rmode, status;
  1632.  
  1633.   // UF occurs
  1634.   if (expon + MAX_FORMAT_DIGITS_128 < 0) {
  1635. #ifdef SET_STATUS_FLAGS
  1636.     __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  1637. #endif
  1638.     pres->w[1] = sgn;
  1639.     pres->w[0] = 0;
  1640. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1641. #ifndef IEEE_ROUND_NEAREST
  1642.     if ((sgn && *prounding_mode == ROUNDING_DOWN)
  1643.         || (!sgn && *prounding_mode == ROUNDING_UP))
  1644.       pres->w[0] = 1ull;
  1645. #endif
  1646. #endif
  1647.     return pres;
  1648.   }
  1649.  
  1650.   ed2 = 0 - expon;
  1651.   // add rounding constant to CQ
  1652. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1653. #ifndef IEEE_ROUND_NEAREST
  1654.   rmode = *prounding_mode;
  1655.   if (sgn && (unsigned) (rmode - 1) < 2)
  1656.     rmode = 3 - rmode;
  1657. #else
  1658.   rmode = 0;
  1659. #endif
  1660. #else
  1661.   rmode = 0;
  1662. #endif
  1663.  
  1664.   T128 = round_const_table_128[rmode][ed2];
  1665.   __add_carry_out (CQ.w[0], carry, T128.w[0], CQ.w[0]);
  1666.   CQ.w[1] = CQ.w[1] + T128.w[1] + carry;
  1667.  
  1668.   TP128 = reciprocals10_128[ed2];
  1669.   __mul_128x128_full (Qh, Ql, CQ, TP128);
  1670.   amount = recip_scale[ed2];
  1671.  
  1672.   if (amount >= 64) {
  1673.     CQ.w[0] = Qh.w[1] >> (amount - 64);
  1674.     CQ.w[1] = 0;
  1675.   } else {
  1676.     __shr_128 (CQ, Qh, amount);
  1677.   }
  1678.  
  1679.   expon = 0;
  1680.  
  1681. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1682. #ifndef IEEE_ROUND_NEAREST
  1683.   if (!(*prounding_mode))
  1684. #endif
  1685.     if (CQ.w[0] & 1) {
  1686.       // check whether fractional part of initial_P/10^ed1 is exactly .5
  1687.  
  1688.       // get remainder
  1689.       __shl_128_long (Qh1, Qh, (128 - amount));
  1690.  
  1691.       if (!Qh1.w[1] && !Qh1.w[0]
  1692.           && (Ql.w[1] < reciprocals10_128[ed2].w[1]
  1693.               || (Ql.w[1] == reciprocals10_128[ed2].w[1]
  1694.                   && Ql.w[0] < reciprocals10_128[ed2].w[0]))) {
  1695.         CQ.w[0]--;
  1696.       }
  1697.     }
  1698. #endif
  1699.  
  1700. #ifdef SET_STATUS_FLAGS
  1701.  
  1702.   if (is_inexact (fpsc))
  1703.     __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
  1704.   else {
  1705.     status = INEXACT_EXCEPTION;
  1706.     // get remainder
  1707.     __shl_128_long (Qh1, Qh, (128 - amount));
  1708.  
  1709.     switch (rmode) {
  1710.     case ROUNDING_TO_NEAREST:
  1711.     case ROUNDING_TIES_AWAY:
  1712.       // test whether fractional part is 0
  1713.       if (Qh1.w[1] == 0x8000000000000000ull && (!Qh1.w[0])
  1714.           && (Ql.w[1] < reciprocals10_128[ed2].w[1]
  1715.               || (Ql.w[1] == reciprocals10_128[ed2].w[1]
  1716.                   && Ql.w[0] < reciprocals10_128[ed2].w[0])))
  1717.         status = EXACT_STATUS;
  1718.       break;
  1719.     case ROUNDING_DOWN:
  1720.     case ROUNDING_TO_ZERO:
  1721.       if ((!Qh1.w[1]) && (!Qh1.w[0])
  1722.           && (Ql.w[1] < reciprocals10_128[ed2].w[1]
  1723.               || (Ql.w[1] == reciprocals10_128[ed2].w[1]
  1724.                   && Ql.w[0] < reciprocals10_128[ed2].w[0])))
  1725.         status = EXACT_STATUS;
  1726.       break;
  1727.     default:
  1728.       // round up
  1729.       __add_carry_out (Stemp.w[0], CY, Ql.w[0],
  1730.                        reciprocals10_128[ed2].w[0]);
  1731.       __add_carry_in_out (Stemp.w[1], carry, Ql.w[1],
  1732.                           reciprocals10_128[ed2].w[1], CY);
  1733.       __shr_128_long (Qh, Qh1, (128 - amount));
  1734.       Tmp.w[0] = 1;
  1735.       Tmp.w[1] = 0;
  1736.       __shl_128_long (Tmp1, Tmp, amount);
  1737.       Qh.w[0] += carry;
  1738.       if (Qh.w[0] < carry)
  1739.         Qh.w[1]++;
  1740.       if (__unsigned_compare_ge_128 (Qh, Tmp1))
  1741.         status = EXACT_STATUS;
  1742.     }
  1743.  
  1744.     if (status != EXACT_STATUS)
  1745.       __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status);
  1746.   }
  1747.  
  1748. #endif
  1749.  
  1750.   pres->w[1] = sgn | CQ.w[1];
  1751.   pres->w[0] = CQ.w[0];
  1752.  
  1753.   return pres;
  1754.  
  1755. }
  1756.  
  1757.  
  1758.  
  1759. //
  1760. //  BID128 unpack, input passed by value
  1761. //
  1762. __BID_INLINE__ UINT64
  1763. unpack_BID128_value (UINT64 * psign_x, int *pexponent_x,
  1764.                      UINT128 * pcoefficient_x, UINT128 x) {
  1765.   UINT128 coeff, T33, T34;
  1766.   UINT64 ex;
  1767.  
  1768.   *psign_x = (x.w[1]) & 0x8000000000000000ull;
  1769.  
  1770.   // special encodings
  1771.   if ((x.w[1] & INFINITY_MASK64) >= SPECIAL_ENCODING_MASK64) {
  1772.     if ((x.w[1] & INFINITY_MASK64) < INFINITY_MASK64) {
  1773.       // non-canonical input
  1774.       pcoefficient_x->w[0] = 0;
  1775.       pcoefficient_x->w[1] = 0;
  1776.       ex = (x.w[1]) >> 47;
  1777.       *pexponent_x = ((int) ex) & EXPONENT_MASK128;
  1778.       return 0;
  1779.     }
  1780.     // 10^33
  1781.     T33 = power10_table_128[33];
  1782.     /*coeff.w[0] = x.w[0];
  1783.        coeff.w[1] = (x.w[1]) & LARGE_COEFF_MASK128;
  1784.        pcoefficient_x->w[0] = x.w[0];
  1785.        pcoefficient_x->w[1] = x.w[1];
  1786.        if (__unsigned_compare_ge_128 (coeff, T33)) // non-canonical
  1787.        pcoefficient_x->w[1] &= (~LARGE_COEFF_MASK128); */
  1788.  
  1789.     pcoefficient_x->w[0] = x.w[0];
  1790.     pcoefficient_x->w[1] = (x.w[1]) & 0x00003fffffffffffull;
  1791.     if (__unsigned_compare_ge_128 ((*pcoefficient_x), T33))     // non-canonical
  1792.     {
  1793.       pcoefficient_x->w[1] = (x.w[1]) & 0xfe00000000000000ull;
  1794.       pcoefficient_x->w[0] = 0;
  1795.     } else
  1796.       pcoefficient_x->w[1] = (x.w[1]) & 0xfe003fffffffffffull;
  1797.     if ((x.w[1] & NAN_MASK64) == INFINITY_MASK64) {
  1798.       pcoefficient_x->w[0] = 0;
  1799.       pcoefficient_x->w[1] = x.w[1] & SINFINITY_MASK64;
  1800.     }
  1801.     *pexponent_x = 0;
  1802.     return 0;   // NaN or Infinity
  1803.   }
  1804.  
  1805.   coeff.w[0] = x.w[0];
  1806.   coeff.w[1] = (x.w[1]) & SMALL_COEFF_MASK128;
  1807.  
  1808.   // 10^34
  1809.   T34 = power10_table_128[34];
  1810.   // check for non-canonical values
  1811.   if (__unsigned_compare_ge_128 (coeff, T34))
  1812.     coeff.w[0] = coeff.w[1] = 0;
  1813.  
  1814.   pcoefficient_x->w[0] = coeff.w[0];
  1815.   pcoefficient_x->w[1] = coeff.w[1];
  1816.  
  1817.   ex = (x.w[1]) >> 49;
  1818.   *pexponent_x = ((int) ex) & EXPONENT_MASK128;
  1819.  
  1820.   return coeff.w[0] | coeff.w[1];
  1821. }
  1822.  
  1823.  
  1824. //
  1825. //  BID128 unpack, input pased by reference
  1826. //
  1827. __BID_INLINE__ UINT64
  1828. unpack_BID128 (UINT64 * psign_x, int *pexponent_x,
  1829.                UINT128 * pcoefficient_x, UINT128 * px) {
  1830.   UINT128 coeff, T33, T34;
  1831.   UINT64 ex;
  1832.  
  1833.   *psign_x = (px->w[1]) & 0x8000000000000000ull;
  1834.  
  1835.   // special encodings
  1836.   if ((px->w[1] & INFINITY_MASK64) >= SPECIAL_ENCODING_MASK64) {
  1837.     if ((px->w[1] & INFINITY_MASK64) < INFINITY_MASK64) {
  1838.       // non-canonical input
  1839.       pcoefficient_x->w[0] = 0;
  1840.       pcoefficient_x->w[1] = 0;
  1841.       ex = (px->w[1]) >> 47;
  1842.       *pexponent_x = ((int) ex) & EXPONENT_MASK128;
  1843.       return 0;
  1844.     }
  1845.     // 10^33
  1846.     T33 = power10_table_128[33];
  1847.     coeff.w[0] = px->w[0];
  1848.     coeff.w[1] = (px->w[1]) & LARGE_COEFF_MASK128;
  1849.     pcoefficient_x->w[0] = px->w[0];
  1850.     pcoefficient_x->w[1] = px->w[1];
  1851.     if (__unsigned_compare_ge_128 (coeff, T33)) {       // non-canonical
  1852.       pcoefficient_x->w[1] &= (~LARGE_COEFF_MASK128);
  1853.       pcoefficient_x->w[0] = 0;
  1854.     }
  1855.     *pexponent_x = 0;
  1856.     return 0;   // NaN or Infinity
  1857.   }
  1858.  
  1859.   coeff.w[0] = px->w[0];
  1860.   coeff.w[1] = (px->w[1]) & SMALL_COEFF_MASK128;
  1861.  
  1862.   // 10^34
  1863.   T34 = power10_table_128[34];
  1864.   // check for non-canonical values
  1865.   if (__unsigned_compare_ge_128 (coeff, T34))
  1866.     coeff.w[0] = coeff.w[1] = 0;
  1867.  
  1868.   pcoefficient_x->w[0] = coeff.w[0];
  1869.   pcoefficient_x->w[1] = coeff.w[1];
  1870.  
  1871.   ex = (px->w[1]) >> 49;
  1872.   *pexponent_x = ((int) ex) & EXPONENT_MASK128;
  1873.  
  1874.   return coeff.w[0] | coeff.w[1];
  1875. }
  1876.  
  1877. //
  1878. //   Pack macro checks for overflow, but not underflow
  1879. //
  1880. __BID_INLINE__ UINT128 *
  1881. get_BID128_very_fast_OF (UINT128 * pres, UINT64 sgn, int expon,
  1882.                          UINT128 coeff, unsigned *prounding_mode,
  1883.                          unsigned *fpsc) {
  1884.   UINT128 T;
  1885.   UINT64 tmp, tmp2;
  1886.  
  1887.   if ((unsigned) expon > DECIMAL_MAX_EXPON_128) {
  1888.  
  1889.     if (expon - MAX_FORMAT_DIGITS_128 <= DECIMAL_MAX_EXPON_128) {
  1890.       T = power10_table_128[MAX_FORMAT_DIGITS_128 - 1];
  1891.       while (__unsigned_compare_gt_128 (T, coeff)
  1892.              && expon > DECIMAL_MAX_EXPON_128) {
  1893.         coeff.w[1] =
  1894.           (coeff.w[1] << 3) + (coeff.w[1] << 1) + (coeff.w[0] >> 61) +
  1895.           (coeff.w[0] >> 63);
  1896.         tmp2 = coeff.w[0] << 3;
  1897.         coeff.w[0] = (coeff.w[0] << 1) + tmp2;
  1898.         if (coeff.w[0] < tmp2)
  1899.           coeff.w[1]++;
  1900.  
  1901.         expon--;
  1902.       }
  1903.     }
  1904.     if ((unsigned) expon > DECIMAL_MAX_EXPON_128) {
  1905.       // OF
  1906. #ifdef SET_STATUS_FLAGS
  1907.       __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  1908. #endif
  1909. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  1910. #ifndef IEEE_ROUND_NEAREST
  1911.       if (*prounding_mode == ROUNDING_TO_ZERO
  1912.           || (sgn && *prounding_mode == ROUNDING_UP) || (!sgn
  1913.                                                          &&
  1914.                                                          *prounding_mode
  1915.                                                          ==
  1916.                                                          ROUNDING_DOWN))
  1917.       {
  1918.         pres->w[1] = sgn | LARGEST_BID128_HIGH;
  1919.         pres->w[0] = LARGEST_BID128_LOW;
  1920.       } else
  1921. #endif
  1922. #endif
  1923.       {
  1924.         pres->w[1] = sgn | INFINITY_MASK64;
  1925.         pres->w[0] = 0;
  1926.       }
  1927.       return pres;
  1928.     }
  1929.   }
  1930.  
  1931.   pres->w[0] = coeff.w[0];
  1932.   tmp = expon;
  1933.   tmp <<= 49;
  1934.   pres->w[1] = sgn | tmp | coeff.w[1];
  1935.  
  1936.   return pres;
  1937. }
  1938.  
  1939.  
  1940. //
  1941. //   No overflow/underflow checks
  1942. //   No checking for coefficient == 10^34 (rounding artifact)
  1943. //
  1944. __BID_INLINE__ UINT128 *
  1945. get_BID128_very_fast (UINT128 * pres, UINT64 sgn, int expon,
  1946.                       UINT128 coeff) {
  1947.   UINT64 tmp;
  1948.  
  1949.   pres->w[0] = coeff.w[0];
  1950.   tmp = expon;
  1951.   tmp <<= 49;
  1952.   pres->w[1] = sgn | tmp | coeff.w[1];
  1953.  
  1954.   return pres;
  1955. }
  1956.  
  1957. //
  1958. //   No overflow/underflow checks
  1959. //
  1960. __BID_INLINE__ UINT128 *
  1961. get_BID128_fast (UINT128 * pres, UINT64 sgn, int expon, UINT128 coeff) {
  1962.   UINT64 tmp;
  1963.  
  1964.   // coeff==10^34?
  1965.   if (coeff.w[1] == 0x0001ed09bead87c0ull
  1966.       && coeff.w[0] == 0x378d8e6400000000ull) {
  1967.     expon++;
  1968.     // set coefficient to 10^33
  1969.     coeff.w[1] = 0x0000314dc6448d93ull;
  1970.     coeff.w[0] = 0x38c15b0a00000000ull;
  1971.   }
  1972.  
  1973.   pres->w[0] = coeff.w[0];
  1974.   tmp = expon;
  1975.   tmp <<= 49;
  1976.   pres->w[1] = sgn | tmp | coeff.w[1];
  1977.  
  1978.   return pres;
  1979. }
  1980.  
  1981. //
  1982. //   General BID128 pack macro
  1983. //
  1984. __BID_INLINE__ UINT128 *
  1985. get_BID128 (UINT128 * pres, UINT64 sgn, int expon, UINT128 coeff,
  1986.             unsigned *prounding_mode, unsigned *fpsc) {
  1987.   UINT128 T;
  1988.   UINT64 tmp, tmp2;
  1989.  
  1990.   // coeff==10^34?
  1991.   if (coeff.w[1] == 0x0001ed09bead87c0ull
  1992.       && coeff.w[0] == 0x378d8e6400000000ull) {
  1993.     expon++;
  1994.     // set coefficient to 10^33
  1995.     coeff.w[1] = 0x0000314dc6448d93ull;
  1996.     coeff.w[0] = 0x38c15b0a00000000ull;
  1997.   }
  1998.   // check OF, UF
  1999.   if (expon < 0 || expon > DECIMAL_MAX_EXPON_128) {
  2000.     // check UF
  2001.     if (expon < 0) {
  2002.       return handle_UF_128 (pres, sgn, expon, coeff, prounding_mode,
  2003.                             fpsc);
  2004.     }
  2005.  
  2006.     if (expon - MAX_FORMAT_DIGITS_128 <= DECIMAL_MAX_EXPON_128) {
  2007.       T = power10_table_128[MAX_FORMAT_DIGITS_128 - 1];
  2008.       while (__unsigned_compare_gt_128 (T, coeff)
  2009.              && expon > DECIMAL_MAX_EXPON_128) {
  2010.         coeff.w[1] =
  2011.           (coeff.w[1] << 3) + (coeff.w[1] << 1) + (coeff.w[0] >> 61) +
  2012.           (coeff.w[0] >> 63);
  2013.         tmp2 = coeff.w[0] << 3;
  2014.         coeff.w[0] = (coeff.w[0] << 1) + tmp2;
  2015.         if (coeff.w[0] < tmp2)
  2016.           coeff.w[1]++;
  2017.  
  2018.         expon--;
  2019.       }
  2020.     }
  2021.     if (expon > DECIMAL_MAX_EXPON_128) {
  2022.       if (!(coeff.w[1] | coeff.w[0])) {
  2023.         pres->w[1] = sgn | (((UINT64) DECIMAL_MAX_EXPON_128) << 49);
  2024.         pres->w[0] = 0;
  2025.         return pres;
  2026.       }
  2027.       // OF
  2028. #ifdef SET_STATUS_FLAGS
  2029.       __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  2030. #endif
  2031. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  2032. #ifndef IEEE_ROUND_NEAREST
  2033.       if (*prounding_mode == ROUNDING_TO_ZERO
  2034.           || (sgn && *prounding_mode == ROUNDING_UP) || (!sgn
  2035.                                                          &&
  2036.                                                          *prounding_mode
  2037.                                                          ==
  2038.                                                          ROUNDING_DOWN))
  2039.       {
  2040.         pres->w[1] = sgn | LARGEST_BID128_HIGH;
  2041.         pres->w[0] = LARGEST_BID128_LOW;
  2042.       } else
  2043. #endif
  2044. #endif
  2045.       {
  2046.         pres->w[1] = sgn | INFINITY_MASK64;
  2047.         pres->w[0] = 0;
  2048.       }
  2049.       return pres;
  2050.     }
  2051.   }
  2052.  
  2053.   pres->w[0] = coeff.w[0];
  2054.   tmp = expon;
  2055.   tmp <<= 49;
  2056.   pres->w[1] = sgn | tmp | coeff.w[1];
  2057.  
  2058.   return pres;
  2059. }
  2060.  
  2061.  
  2062. //
  2063. //  Macro used for conversions from string
  2064. //        (no additional arguments given for rounding mode, status flags)
  2065. //
  2066. __BID_INLINE__ UINT128 *
  2067. get_BID128_string (UINT128 * pres, UINT64 sgn, int expon, UINT128 coeff) {
  2068.   UINT128 D2, D8;
  2069.   UINT64 tmp;
  2070.   unsigned rmode = 0, status;
  2071.  
  2072.   // coeff==10^34?
  2073.   if (coeff.w[1] == 0x0001ed09bead87c0ull
  2074.       && coeff.w[0] == 0x378d8e6400000000ull) {
  2075.     expon++;
  2076.     // set coefficient to 10^33
  2077.     coeff.w[1] = 0x0000314dc6448d93ull;
  2078.     coeff.w[0] = 0x38c15b0a00000000ull;
  2079.   }
  2080.   // check OF, UF
  2081.   if ((unsigned) expon > DECIMAL_MAX_EXPON_128) {
  2082.     // check UF
  2083.     if (expon < 0)
  2084.       return handle_UF_128 (pres, sgn, expon, coeff, &rmode, &status);
  2085.  
  2086.     // OF
  2087.  
  2088.     if (expon < DECIMAL_MAX_EXPON_128 + 34) {
  2089.       while (expon > DECIMAL_MAX_EXPON_128 &&
  2090.              (coeff.w[1] < power10_table_128[33].w[1] ||
  2091.               (coeff.w[1] == power10_table_128[33].w[1]
  2092.                && coeff.w[0] < power10_table_128[33].w[0]))) {
  2093.         D2.w[1] = (coeff.w[1] << 1) | (coeff.w[0] >> 63);
  2094.         D2.w[0] = coeff.w[0] << 1;
  2095.         D8.w[1] = (coeff.w[1] << 3) | (coeff.w[0] >> 61);
  2096.         D8.w[0] = coeff.w[0] << 3;
  2097.  
  2098.         __add_128_128 (coeff, D2, D8);
  2099.         expon--;
  2100.       }
  2101.     } else if (!(coeff.w[0] | coeff.w[1]))
  2102.       expon = DECIMAL_MAX_EXPON_128;
  2103.  
  2104.     if (expon > DECIMAL_MAX_EXPON_128) {
  2105.       pres->w[1] = sgn | INFINITY_MASK64;
  2106.       pres->w[0] = 0;
  2107.       switch (rmode) {
  2108.       case ROUNDING_DOWN:
  2109.         if (!sgn) {
  2110.           pres->w[1] = LARGEST_BID128_HIGH;
  2111.           pres->w[0] = LARGEST_BID128_LOW;
  2112.         }
  2113.         break;
  2114.       case ROUNDING_TO_ZERO:
  2115.         pres->w[1] = sgn | LARGEST_BID128_HIGH;
  2116.         pres->w[0] = LARGEST_BID128_LOW;
  2117.         break;
  2118.       case ROUNDING_UP:
  2119.         // round up
  2120.         if (sgn) {
  2121.           pres->w[1] = sgn | LARGEST_BID128_HIGH;
  2122.           pres->w[0] = LARGEST_BID128_LOW;
  2123.         }
  2124.         break;
  2125.       }
  2126.  
  2127.       return pres;
  2128.     }
  2129.   }
  2130.  
  2131.   pres->w[0] = coeff.w[0];
  2132.   tmp = expon;
  2133.   tmp <<= 49;
  2134.   pres->w[1] = sgn | tmp | coeff.w[1];
  2135.  
  2136.   return pres;
  2137. }
  2138.  
  2139.  
  2140.  
  2141. /*****************************************************************************
  2142. *
  2143. *    BID32 pack/unpack macros
  2144. *
  2145. *****************************************************************************/
  2146.  
  2147.  
  2148. __BID_INLINE__ UINT32
  2149. unpack_BID32 (UINT32 * psign_x, int *pexponent_x,
  2150.               UINT32 * pcoefficient_x, UINT32 x) {
  2151.   UINT32 tmp;
  2152.  
  2153.   *psign_x = x & 0x80000000;
  2154.  
  2155.   if ((x & SPECIAL_ENCODING_MASK32) == SPECIAL_ENCODING_MASK32) {
  2156.     // special encodings
  2157.     if ((x & INFINITY_MASK32) == INFINITY_MASK32) {
  2158.       *pcoefficient_x = x & 0xfe0fffff;
  2159.       if ((x & 0x000fffff) >= 1000000)
  2160.         *pcoefficient_x = x & 0xfe000000;
  2161.       if ((x & NAN_MASK32) == INFINITY_MASK32)
  2162.         *pcoefficient_x = x & 0xf8000000;
  2163.       *pexponent_x = 0;
  2164.       return 0; // NaN or Infinity
  2165.     }
  2166.     // coefficient
  2167.     *pcoefficient_x = (x & SMALL_COEFF_MASK32) | LARGE_COEFF_HIGH_BIT32;
  2168.     // check for non-canonical value
  2169.     if (*pcoefficient_x >= 10000000)
  2170.       *pcoefficient_x = 0;
  2171.     // get exponent
  2172.     tmp = x >> 21;
  2173.     *pexponent_x = tmp & EXPONENT_MASK32;
  2174.     return 1;
  2175.   }
  2176.   // exponent
  2177.   tmp = x >> 23;
  2178.   *pexponent_x = tmp & EXPONENT_MASK32;
  2179.   // coefficient
  2180.   *pcoefficient_x = (x & LARGE_COEFF_MASK32);
  2181.  
  2182.   return *pcoefficient_x;
  2183. }
  2184.  
  2185. //
  2186. //   General pack macro for BID32
  2187. //
  2188. __BID_INLINE__ UINT32
  2189. get_BID32 (UINT32 sgn, int expon, UINT64 coeff, int rmode,
  2190.            unsigned *fpsc) {
  2191.   UINT128 Q;
  2192.   UINT64 C64, remainder_h, carry, Stemp;
  2193.   UINT32 r, mask;
  2194.   int extra_digits, amount, amount2;
  2195.   unsigned status;
  2196.  
  2197.   if (coeff > 9999999ull) {
  2198.     expon++;
  2199.     coeff = 1000000ull;
  2200.   }
  2201.   // check for possible underflow/overflow
  2202.   if (((unsigned) expon) > DECIMAL_MAX_EXPON_32) {
  2203.     if (expon < 0) {
  2204.       // underflow
  2205.       if (expon + MAX_FORMAT_DIGITS_32 < 0) {
  2206. #ifdef SET_STATUS_FLAGS
  2207.         __set_status_flags (fpsc,
  2208.                             UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  2209. #endif
  2210. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  2211. #ifndef IEEE_ROUND_NEAREST
  2212.         if (rmode == ROUNDING_DOWN && sgn)
  2213.           return 0x80000001;
  2214.         if (rmode == ROUNDING_UP && !sgn)
  2215.           return 1;
  2216. #endif
  2217. #endif
  2218.         // result is 0
  2219.         return sgn;
  2220.       }
  2221.       // get digits to be shifted out
  2222. #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
  2223.       rmode = 0;
  2224. #endif
  2225. #ifdef IEEE_ROUND_NEAREST
  2226.       rmode = 0;
  2227. #endif
  2228. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  2229. #ifndef IEEE_ROUND_NEAREST
  2230.       if (sgn && (unsigned) (rmode - 1) < 2)
  2231.         rmode = 3 - rmode;
  2232. #endif
  2233. #endif
  2234.  
  2235.       extra_digits = -expon;
  2236.       coeff += round_const_table[rmode][extra_digits];
  2237.  
  2238.       // get coeff*(2^M[extra_digits])/10^extra_digits
  2239.       __mul_64x64_to_128 (Q, coeff, reciprocals10_64[extra_digits]);
  2240.  
  2241.       // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
  2242.       amount = short_recip_scale[extra_digits];
  2243.  
  2244.       C64 = Q.w[1] >> amount;
  2245.  
  2246. #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
  2247. #ifndef IEEE_ROUND_NEAREST
  2248.       if (rmode == 0)   //ROUNDING_TO_NEAREST
  2249. #endif
  2250.         if (C64 & 1) {
  2251.           // check whether fractional part of initial_P/10^extra_digits is exactly .5
  2252.  
  2253.           // get remainder
  2254.           amount2 = 64 - amount;
  2255.           remainder_h = 0;
  2256.           remainder_h--;
  2257.           remainder_h >>= amount2;
  2258.           remainder_h = remainder_h & Q.w[1];
  2259.  
  2260.           if (!remainder_h && (Q.w[0] < reciprocals10_64[extra_digits])) {
  2261.             C64--;
  2262.           }
  2263.         }
  2264. #endif
  2265.  
  2266. #ifdef SET_STATUS_FLAGS
  2267.  
  2268.       if (is_inexact (fpsc))
  2269.         __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
  2270.       else {
  2271.         status = INEXACT_EXCEPTION;
  2272.         // get remainder
  2273.         remainder_h = Q.w[1] << (64 - amount);
  2274.  
  2275.         switch (rmode) {
  2276.         case ROUNDING_TO_NEAREST:
  2277.         case ROUNDING_TIES_AWAY:
  2278.           // test whether fractional part is 0
  2279.           if (remainder_h == 0x8000000000000000ull
  2280.               && (Q.w[0] < reciprocals10_64[extra_digits]))
  2281.             status = EXACT_STATUS;
  2282.           break;
  2283.         case ROUNDING_DOWN:
  2284.         case ROUNDING_TO_ZERO:
  2285.           if (!remainder_h && (Q.w[0] < reciprocals10_64[extra_digits]))
  2286.             status = EXACT_STATUS;
  2287.           break;
  2288.         default:
  2289.           // round up
  2290.           __add_carry_out (Stemp, carry, Q.w[0],
  2291.                            reciprocals10_64[extra_digits]);
  2292.           if ((remainder_h >> (64 - amount)) + carry >=
  2293.               (((UINT64) 1) << amount))
  2294.             status = EXACT_STATUS;
  2295.         }
  2296.  
  2297.         if (status != EXACT_STATUS)
  2298.           __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status);
  2299.       }
  2300.  
  2301. #endif
  2302.  
  2303.       return sgn | (UINT32) C64;
  2304.     }
  2305.  
  2306.     while (coeff < 1000000 && expon > DECIMAL_MAX_EXPON_32) {
  2307.       coeff = (coeff << 3) + (coeff << 1);
  2308.       expon--;
  2309.     }
  2310.     if (((unsigned) expon) > DECIMAL_MAX_EXPON_32) {
  2311.       __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
  2312.       // overflow
  2313.       r = sgn | INFINITY_MASK32;
  2314.       switch (rmode) {
  2315.       case ROUNDING_DOWN:
  2316.         if (!sgn)
  2317.           r = LARGEST_BID32;
  2318.         break;
  2319.       case ROUNDING_TO_ZERO:
  2320.         r = sgn | LARGEST_BID32;
  2321.         break;
  2322.       case ROUNDING_UP:
  2323.         // round up
  2324.         if (sgn)
  2325.           r = sgn | LARGEST_BID32;
  2326.       }
  2327.       return r;
  2328.     }
  2329.   }
  2330.  
  2331.   mask = 1 << 23;
  2332.  
  2333.   // check whether coefficient fits in DECIMAL_COEFF_FIT bits
  2334.   if (coeff < mask) {
  2335.     r = expon;
  2336.     r <<= 23;
  2337.     r |= ((UINT32) coeff | sgn);
  2338.     return r;
  2339.   }
  2340.   // special format
  2341.  
  2342.   r = expon;
  2343.   r <<= 21;
  2344.   r |= (sgn | SPECIAL_ENCODING_MASK32);
  2345.   // add coeff, without leading bits
  2346.   mask = (1 << 21) - 1;
  2347.   r |= (((UINT32) coeff) & mask);
  2348.  
  2349.   return r;
  2350. }
  2351.  
  2352.  
  2353.  
  2354. //
  2355. //   no overflow/underflow checks
  2356. //
  2357. __BID_INLINE__ UINT32
  2358. very_fast_get_BID32 (UINT32 sgn, int expon, UINT32 coeff) {
  2359.   UINT32 r, mask;
  2360.  
  2361.   mask = 1 << 23;
  2362.  
  2363.   // check whether coefficient fits in 10*2+3 bits
  2364.   if (coeff < mask) {
  2365.     r = expon;
  2366.     r <<= 23;
  2367.     r |= (coeff | sgn);
  2368.     return r;
  2369.   }
  2370.   // special format
  2371.   r = expon;
  2372.   r <<= 21;
  2373.   r |= (sgn | SPECIAL_ENCODING_MASK32);
  2374.   // add coeff, without leading bits
  2375.   mask = (1 << 21) - 1;
  2376.   coeff &= mask;
  2377.   r |= coeff;
  2378.  
  2379.   return r;
  2380. }
  2381.  
  2382.  
  2383.  
  2384. /*************************************************************
  2385.  *
  2386.  *************************************************************/
  2387. typedef
  2388. ALIGN (16)
  2389.      struct {
  2390.        UINT64 w[6];
  2391.      } UINT384;
  2392.      typedef ALIGN (16)
  2393.      struct {
  2394.        UINT64 w[8];
  2395.      } UINT512;
  2396.  
  2397. // #define P                               34
  2398. #define MASK_STEERING_BITS              0x6000000000000000ull
  2399. #define MASK_BINARY_EXPONENT1           0x7fe0000000000000ull
  2400. #define MASK_BINARY_SIG1                0x001fffffffffffffull
  2401. #define MASK_BINARY_EXPONENT2           0x1ff8000000000000ull
  2402.     //used to take G[2:w+3] (sec 3.3)
  2403. #define MASK_BINARY_SIG2                0x0007ffffffffffffull
  2404.     //used to mask out G4:T0 (sec 3.3)
  2405. #define MASK_BINARY_OR2                 0x0020000000000000ull
  2406.     //used to prefix 8+G4 to T (sec 3.3)
  2407. #define UPPER_EXPON_LIMIT               51
  2408. #define MASK_EXP                        0x7ffe000000000000ull
  2409. #define MASK_SPECIAL                    0x7800000000000000ull
  2410. #define MASK_NAN                        0x7c00000000000000ull
  2411. #define MASK_SNAN                       0x7e00000000000000ull
  2412. #define MASK_ANY_INF                    0x7c00000000000000ull
  2413. #define MASK_INF                        0x7800000000000000ull
  2414. #define MASK_SIGN                       0x8000000000000000ull
  2415. #define MASK_COEFF                      0x0001ffffffffffffull
  2416. #define BIN_EXP_BIAS                    (0x1820ull << 49)
  2417.  
  2418. #define EXP_MIN                         0x0000000000000000ull
  2419.    // EXP_MIN = (-6176 + 6176) << 49
  2420. #define EXP_MAX                         0x5ffe000000000000ull
  2421.   // EXP_MAX = (6111 + 6176) << 49
  2422. #define EXP_MAX_P1                      0x6000000000000000ull
  2423.   // EXP_MAX + 1 = (6111 + 6176 + 1) << 49
  2424. #define EXP_P1                          0x0002000000000000ull
  2425.   // EXP_ P1= 1 << 49
  2426. #define expmin                            -6176
  2427.   // min unbiased exponent
  2428. #define expmax                            6111
  2429.   // max unbiased exponent
  2430. #define expmin16                          -398
  2431.   // min unbiased exponent
  2432. #define expmax16                          369
  2433.   // max unbiased exponent
  2434.  
  2435. #define SIGNMASK32 0x80000000
  2436. #define BID64_SIG_MAX 0x002386F26FC0ffffull
  2437. #define SIGNMASK64    0x8000000000000000ull
  2438.  
  2439. // typedef unsigned int FPSC; // floating-point status and control
  2440.         // bit31:
  2441.         // bit30:
  2442.         // bit29:
  2443.         // bit28:
  2444.         // bit27:
  2445.         // bit26:
  2446.         // bit25:
  2447.         // bit24:
  2448.         // bit23:
  2449.         // bit22:
  2450.         // bit21:
  2451.         // bit20:
  2452.         // bit19:
  2453.         // bit18:
  2454.         // bit17:
  2455.         // bit16:
  2456.         // bit15:
  2457.         // bit14: RC:2
  2458.         // bit13: RC:1
  2459.         // bit12: RC:0
  2460.         // bit11: PM
  2461.         // bit10: UM
  2462.         // bit9:  OM
  2463.         // bit8:  ZM
  2464.         // bit7:  DM
  2465.         // bit6:  IM
  2466.         // bit5:  PE
  2467.         // bit4:  UE
  2468.         // bit3:  OE
  2469.         // bit2:  ZE
  2470.         // bit1:  DE
  2471.         // bit0:  IE
  2472.  
  2473. #define ROUNDING_MODE_MASK      0x00007000
  2474.  
  2475.      typedef struct _DEC_DIGITS {
  2476.        unsigned int digits;
  2477.        UINT64 threshold_hi;
  2478.        UINT64 threshold_lo;
  2479.        unsigned int digits1;
  2480.      } DEC_DIGITS;
  2481.  
  2482.      extern DEC_DIGITS nr_digits[];
  2483.      extern UINT64 midpoint64[];
  2484.      extern UINT128 midpoint128[];
  2485.      extern UINT192 midpoint192[];
  2486.      extern UINT256 midpoint256[];
  2487.      extern UINT64 ten2k64[];
  2488.      extern UINT128 ten2k128[];
  2489.      extern UINT256 ten2k256[];
  2490.      extern UINT128 ten2mk128[];
  2491.      extern UINT64 ten2mk64[];
  2492.      extern UINT128 ten2mk128trunc[];
  2493.      extern int shiftright128[];
  2494.      extern UINT64 maskhigh128[];
  2495.      extern UINT64 maskhigh128M[];
  2496.      extern UINT64 maskhigh192M[];
  2497.      extern UINT64 maskhigh256M[];
  2498.      extern UINT64 onehalf128[];
  2499.      extern UINT64 onehalf128M[];
  2500.      extern UINT64 onehalf192M[];
  2501.      extern UINT64 onehalf256M[];
  2502.      extern UINT128 ten2mk128M[];
  2503.      extern UINT128 ten2mk128truncM[];
  2504.      extern UINT192 ten2mk192truncM[];
  2505.      extern UINT256 ten2mk256truncM[];
  2506.      extern int shiftright128M[];
  2507.      extern int shiftright192M[];
  2508.      extern int shiftright256M[];
  2509.      extern UINT192 ten2mk192M[];
  2510.      extern UINT256 ten2mk256M[];
  2511.      extern unsigned char char_table2[];
  2512.      extern unsigned char char_table3[];
  2513.  
  2514.      extern UINT64 ten2m3k64[];
  2515.      extern unsigned int shift_ten2m3k64[];
  2516.      extern UINT128 ten2m3k128[];
  2517.      extern unsigned int shift_ten2m3k128[];
  2518.  
  2519.  
  2520.  
  2521. /***************************************************************************
  2522.  *************** TABLES FOR GENERAL ROUNDING FUNCTIONS *********************
  2523.  ***************************************************************************/
  2524.  
  2525.      extern UINT64 Kx64[];
  2526.      extern unsigned int Ex64m64[];
  2527.      extern UINT64 half64[];
  2528.      extern UINT64 mask64[];
  2529.      extern UINT64 ten2mxtrunc64[];
  2530.  
  2531.      extern UINT128 Kx128[];
  2532.      extern unsigned int Ex128m128[];
  2533.      extern UINT64 half128[];
  2534.      extern UINT64 mask128[];
  2535.      extern UINT128 ten2mxtrunc128[];
  2536.  
  2537.      extern UINT192 Kx192[];
  2538.      extern unsigned int Ex192m192[];
  2539.      extern UINT64 half192[];
  2540.      extern UINT64 mask192[];
  2541.      extern UINT192 ten2mxtrunc192[];
  2542.  
  2543.      extern UINT256 Kx256[];
  2544.      extern unsigned int Ex256m256[];
  2545.      extern UINT64 half256[];
  2546.      extern UINT64 mask256[];
  2547.      extern UINT256 ten2mxtrunc256[];
  2548.  
  2549.      typedef union __bid64_128 {
  2550.        UINT64 b64;
  2551.        UINT128 b128;
  2552.      } BID64_128;
  2553.  
  2554.      BID64_128 bid_fma (unsigned int P0,
  2555.                         BID64_128 x1, unsigned int P1,
  2556.                         BID64_128 y1, unsigned int P2,
  2557.                         BID64_128 z1, unsigned int P3,
  2558.                         unsigned int rnd_mode, FPSC * fpsc);
  2559.  
  2560. #define         P16     16
  2561. #define         P34     34
  2562.  
  2563.      union __int_double {
  2564.        UINT64 i;
  2565.        double d;
  2566.      };
  2567.      typedef union __int_double int_double;
  2568.  
  2569.  
  2570.      union __int_float {
  2571.        UINT32 i;
  2572.        float d;
  2573.      };
  2574.      typedef union __int_float int_float;
  2575.  
  2576. #define SWAP(A,B,T) {\
  2577.         T = A; \
  2578.         A = B; \
  2579.         B = T; \
  2580. }
  2581.  
  2582. // this macro will find coefficient_x to be in [2^A, 2^(A+1) )
  2583. // ie it knows that it is A bits long
  2584. #define NUMBITS(A, coefficient_x, tempx){\
  2585.       temp_x.d=(float)coefficient_x;\
  2586.       A=((tempx.i >>23) & EXPONENT_MASK32) - 0x7f;\
  2587. }
  2588.  
  2589.      enum class_types {
  2590.        signalingNaN,
  2591.        quietNaN,
  2592.        negativeInfinity,
  2593.        negativeNormal,
  2594.        negativeSubnormal,
  2595.        negativeZero,
  2596.        positiveZero,
  2597.        positiveSubnormal,
  2598.        positiveNormal,
  2599.        positiveInfinity
  2600.      };
  2601.  
  2602.      typedef union {
  2603.        UINT64 ui64;
  2604.        double d;
  2605.      } BID_UI64DOUBLE;
  2606.  
  2607. #endif
  2608.