"Fossies" - the Fresh Open Source Software Archive

Member "mpfr-4.0.2/src/get_d64.c" (7 Jan 2019, 19663 Bytes) of package /linux/misc/mpfr-4.0.2.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "get_d64.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 4.0.1_vs_4.0.2.

    1 /* mpfr_get_decimal64 -- convert a multiple precision floating-point number
    2                          to an IEEE 754-2008 decimal64 float
    3 
    4 See https://gcc.gnu.org/ml/gcc/2006-06/msg00691.html,
    5 https://gcc.gnu.org/onlinedocs/gcc/Decimal-Float.html,
    6 and TR 24732 <http://www.open-std.org/jtc1/sc22/wg14/www/projects#24732>.
    7 
    8 Copyright 2006-2019 Free Software Foundation, Inc.
    9 Contributed by the AriC and Caramba projects, INRIA.
   10 
   11 This file is part of the GNU MPFR Library.
   12 
   13 The GNU MPFR Library is free software; you can redistribute it and/or modify
   14 it under the terms of the GNU Lesser General Public License as published by
   15 the Free Software Foundation; either version 3 of the License, or (at your
   16 option) any later version.
   17 
   18 The GNU MPFR Library is distributed in the hope that it will be useful, but
   19 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
   20 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
   21 License for more details.
   22 
   23 You should have received a copy of the GNU Lesser General Public License
   24 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
   25 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
   26 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
   27 
   28 #include "mpfr-impl.h"
   29 
   30 #define ISDIGIT(c) ('0' <= c && c <= '9')
   31 
   32 #ifdef MPFR_WANT_DECIMAL_FLOATS
   33 
   34 #if _MPFR_IEEE_FLOATS
   35 #else
   36 #include "ieee_floats.h"
   37 #endif
   38 
   39 #ifndef DEC64_MAX
   40 # define DEC64_MAX 9.999999999999999E384dd
   41 #endif
   42 
   43 #ifdef DPD_FORMAT
   44 static const int T[1000] = {
   45   0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 32,
   46   33, 34, 35, 36, 37, 38, 39, 40, 41, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57,
   47   64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 80, 81, 82, 83, 84, 85, 86, 87, 88,
   48   89, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 112, 113, 114, 115, 116,
   49   117, 118, 119, 120, 121, 10, 11, 42, 43, 74, 75, 106, 107, 78, 79, 26, 27,
   50   58, 59, 90, 91, 122, 123, 94, 95, 128, 129, 130, 131, 132, 133, 134, 135,
   51   136, 137, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 160, 161, 162,
   52   163, 164, 165, 166, 167, 168, 169, 176, 177, 178, 179, 180, 181, 182, 183,
   53   184, 185, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 208, 209, 210,
   54   211, 212, 213, 214, 215, 216, 217, 224, 225, 226, 227, 228, 229, 230, 231,
   55   232, 233, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 138, 139, 170,
   56   171, 202, 203, 234, 235, 206, 207, 154, 155, 186, 187, 218, 219, 250, 251,
   57   222, 223, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 272, 273, 274,
   58   275, 276, 277, 278, 279, 280, 281, 288, 289, 290, 291, 292, 293, 294, 295,
   59   296, 297, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 320, 321, 322,
   60   323, 324, 325, 326, 327, 328, 329, 336, 337, 338, 339, 340, 341, 342, 343,
   61   344, 345, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 368, 369, 370,
   62   371, 372, 373, 374, 375, 376, 377, 266, 267, 298, 299, 330, 331, 362, 363,
   63   334, 335, 282, 283, 314, 315, 346, 347, 378, 379, 350, 351, 384, 385, 386,
   64   387, 388, 389, 390, 391, 392, 393, 400, 401, 402, 403, 404, 405, 406, 407,
   65   408, 409, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 432, 433, 434,
   66   435, 436, 437, 438, 439, 440, 441, 448, 449, 450, 451, 452, 453, 454, 455,
   67   456, 457, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 480, 481, 482,
   68   483, 484, 485, 486, 487, 488, 489, 496, 497, 498, 499, 500, 501, 502, 503,
   69   504, 505, 394, 395, 426, 427, 458, 459, 490, 491, 462, 463, 410, 411, 442,
   70   443, 474, 475, 506, 507, 478, 479, 512, 513, 514, 515, 516, 517, 518, 519,
   71   520, 521, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 544, 545, 546,
   72   547, 548, 549, 550, 551, 552, 553, 560, 561, 562, 563, 564, 565, 566, 567,
   73   568, 569, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 592, 593, 594,
   74   595, 596, 597, 598, 599, 600, 601, 608, 609, 610, 611, 612, 613, 614, 615,
   75   616, 617, 624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 522, 523, 554,
   76   555, 586, 587, 618, 619, 590, 591, 538, 539, 570, 571, 602, 603, 634, 635,
   77   606, 607, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 656, 657, 658,
   78   659, 660, 661, 662, 663, 664, 665, 672, 673, 674, 675, 676, 677, 678, 679,
   79   680, 681, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 704, 705, 706,
   80   707, 708, 709, 710, 711, 712, 713, 720, 721, 722, 723, 724, 725, 726, 727,
   81   728, 729, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 752, 753, 754,
   82   755, 756, 757, 758, 759, 760, 761, 650, 651, 682, 683, 714, 715, 746, 747,
   83   718, 719, 666, 667, 698, 699, 730, 731, 762, 763, 734, 735, 768, 769, 770,
   84   771, 772, 773, 774, 775, 776, 777, 784, 785, 786, 787, 788, 789, 790, 791,
   85   792, 793, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 816, 817, 818,
   86   819, 820, 821, 822, 823, 824, 825, 832, 833, 834, 835, 836, 837, 838, 839,
   87   840, 841, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 864, 865, 866,
   88   867, 868, 869, 870, 871, 872, 873, 880, 881, 882, 883, 884, 885, 886, 887,
   89   888, 889, 778, 779, 810, 811, 842, 843, 874, 875, 846, 847, 794, 795, 826,
   90   827, 858, 859, 890, 891, 862, 863, 896, 897, 898, 899, 900, 901, 902, 903,
   91   904, 905, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 928, 929, 930,
   92   931, 932, 933, 934, 935, 936, 937, 944, 945, 946, 947, 948, 949, 950, 951,
   93   952, 953, 960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 976, 977, 978,
   94   979, 980, 981, 982, 983, 984, 985, 992, 993, 994, 995, 996, 997, 998, 999,
   95   1000, 1001, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 906,
   96   907, 938, 939, 970, 971, 1002, 1003, 974, 975, 922, 923, 954, 955, 986,
   97   987, 1018, 1019, 990, 991, 12, 13, 268, 269, 524, 525, 780, 781, 46, 47, 28,
   98   29, 284, 285, 540, 541, 796, 797, 62, 63, 44, 45, 300, 301, 556, 557, 812,
   99   813, 302, 303, 60, 61, 316, 317, 572, 573, 828, 829, 318, 319, 76, 77,
  100   332, 333, 588, 589, 844, 845, 558, 559, 92, 93, 348, 349, 604, 605, 860,
  101   861, 574, 575, 108, 109, 364, 365, 620, 621, 876, 877, 814, 815, 124, 125,
  102   380, 381, 636, 637, 892, 893, 830, 831, 14, 15, 270, 271, 526, 527, 782,
  103   783, 110, 111, 30, 31, 286, 287, 542, 543, 798, 799, 126, 127, 140, 141,
  104   396, 397, 652, 653, 908, 909, 174, 175, 156, 157, 412, 413, 668, 669, 924,
  105   925, 190, 191, 172, 173, 428, 429, 684, 685, 940, 941, 430, 431, 188, 189,
  106   444, 445, 700, 701, 956, 957, 446, 447, 204, 205, 460, 461, 716, 717, 972,
  107   973, 686, 687, 220, 221, 476, 477, 732, 733, 988, 989, 702, 703, 236, 237,
  108   492, 493, 748, 749, 1004, 1005, 942, 943, 252, 253, 508, 509, 764, 765,
  109   1020, 1021, 958, 959, 142, 143, 398, 399, 654, 655, 910, 911, 238, 239, 158,
  110   159, 414, 415, 670, 671, 926, 927, 254, 255};
  111 #endif
  112 
  113 /* construct a decimal64 NaN */
  114 static _Decimal64
  115 get_decimal64_nan (void)
  116 {
  117 #if _MPFR_IEEE_FLOATS
  118   union mpfr_ieee_double_extract x;
  119   union ieee_double_decimal64 y;
  120 
  121   x.s.exp = 1984; /* G[0]..G[4] = 11111: quiet NaN */
  122   y.d = x.d;
  123   return y.d64;
  124 #else
  125   return (_Decimal64) MPFR_DBL_NAN;
  126 #endif
  127 }
  128 
  129 /* construct the decimal64 Inf with given sign */
  130 static _Decimal64
  131 get_decimal64_inf (int negative)
  132 {
  133 #if _MPFR_IEEE_FLOATS
  134   union mpfr_ieee_double_extract x;
  135   union ieee_double_decimal64 y;
  136 
  137   x.s.sig = (negative) ? 1 : 0;
  138   x.s.exp = 1920; /* G[0]..G[4] = 11110: Inf */
  139   y.d = x.d;
  140   return y.d64;
  141 #else
  142   return (_Decimal64) (negative ? MPFR_DBL_INFM : MPFR_DBL_INFP);
  143 #endif
  144 }
  145 
  146 /* construct the decimal64 zero with given sign */
  147 static _Decimal64
  148 get_decimal64_zero (int negative)
  149 {
  150   union ieee_double_decimal64 y;
  151 
  152   /* zero has the same representation in binary64 and decimal64 */
  153   y.d = negative ? DBL_NEG_ZERO : 0.0;
  154   return y.d64;
  155 }
  156 
  157 /* construct the decimal64 smallest non-zero with given sign */
  158 static _Decimal64
  159 get_decimal64_min (int negative)
  160 {
  161   return negative ? - 1E-398dd : 1E-398dd;
  162 }
  163 
  164 /* construct the decimal64 largest finite number with given sign */
  165 static _Decimal64
  166 get_decimal64_max (int negative)
  167 {
  168   return negative ? - DEC64_MAX : DEC64_MAX;
  169 }
  170 
  171 /* one-to-one conversion:
  172    s is a decimal string representing a number x = m * 10^e which must be
  173    exactly representable in the decimal64 format, i.e.
  174    (a) the mantissa m has at most 16 decimal digits
  175    (b1) -383 <= e <= 384 with m integer multiple of 10^(-15), |m| < 10
  176    (b2) or -398 <= e <= 369 with m integer, |m| < 10^16.
  177    Assumes s is neither NaN nor +Inf nor -Inf.
  178 */
  179 #if _MPFR_IEEE_FLOATS
  180 static _Decimal64
  181 string_to_Decimal64 (char *s)
  182 {
  183   long int exp = 0;
  184   char m[17];
  185   long n = 0; /* mantissa length */
  186   char *endptr[1];
  187   union mpfr_ieee_double_extract x;
  188   union ieee_double_decimal64 y;
  189 #ifdef DPD_FORMAT
  190   unsigned int G, d1, d2, d3, d4, d5;
  191 #endif
  192 
  193   /* read sign */
  194   if (*s == '-')
  195     {
  196       x.s.sig = 1;
  197       s ++;
  198     }
  199   else
  200     x.s.sig = 0;
  201   /* read mantissa */
  202   while (ISDIGIT (*s))
  203     m[n++] = *s++;
  204   exp = n;
  205   if (*s == '.')
  206     {
  207       s ++;
  208       while (ISDIGIT (*s))
  209         m[n++] = *s++;
  210     }
  211   /* we have exp digits before decimal point, and a total of n digits */
  212   exp -= n; /* we will consider an integer mantissa */
  213   MPFR_ASSERTN(n <= 16);
  214   if (*s == 'E' || *s == 'e')
  215     exp += strtol (s + 1, endptr, 10);
  216   else
  217     *endptr = s;
  218   MPFR_ASSERTN(**endptr == '\0');
  219   MPFR_ASSERTN(-398 <= exp && exp <= (long) (385 - n));
  220   while (n < 16)
  221     {
  222       m[n++] = '0';
  223       exp --;
  224     }
  225   /* now n=16 and -398 <= exp <= 369 */
  226   m[n] = '\0';
  227 
  228   /* compute biased exponent */
  229   exp += 398;
  230 
  231   MPFR_ASSERTN(exp >= -15);
  232   if (exp < 0)
  233     {
  234       int i;
  235       n = -exp;
  236       /* check the last n digits of the mantissa are zero */
  237       for (i = 1; i <= n; i++)
  238         MPFR_ASSERTN(m[16 - n] == '0');
  239       /* shift the first (16-n) digits to the right */
  240       for (i = 16 - n - 1; i >= 0; i--)
  241         m[i + n] = m[i];
  242       /* zero the first n digits */
  243       for (i = 0; i < n; i ++)
  244         m[i] = '0';
  245       exp = 0;
  246     }
  247 
  248   /* now convert to DPD or BID */
  249 #ifdef DPD_FORMAT
  250 #define CH(d) (d - '0')
  251   if (m[0] >= '8')
  252     G = (3 << 11) | ((exp & 768) << 1) | ((CH(m[0]) & 1) << 8);
  253   else
  254     G = ((exp & 768) << 3) | (CH(m[0]) << 8);
  255   /* now the most 5 significant bits of G are filled */
  256   G |= exp & 255;
  257   d1 = T[100 * CH(m[1]) + 10 * CH(m[2]) + CH(m[3])]; /* 10-bit encoding */
  258   d2 = T[100 * CH(m[4]) + 10 * CH(m[5]) + CH(m[6])]; /* 10-bit encoding */
  259   d3 = T[100 * CH(m[7]) + 10 * CH(m[8]) + CH(m[9])]; /* 10-bit encoding */
  260   d4 = T[100 * CH(m[10]) + 10 * CH(m[11]) + CH(m[12])]; /* 10-bit encoding */
  261   d5 = T[100 * CH(m[13]) + 10 * CH(m[14]) + CH(m[15])]; /* 10-bit encoding */
  262   x.s.exp = G >> 2;
  263   x.s.manh = ((G & 3) << 18) | (d1 << 8) | (d2 >> 2);
  264   x.s.manl = (d2 & 3) << 30;
  265   x.s.manl |= (d3 << 20) | (d4 << 10) | d5;
  266 #else /* BID format */
  267   {
  268     mp_size_t rn;
  269     mp_limb_t rp[2];
  270     int case_i = strcmp (m, "9007199254740992") < 0;
  271 
  272     for (n = 0; n < 16; n++)
  273       m[n] -= '0';
  274     rn = mpn_set_str (rp, (unsigned char *) m, 16, 10);
  275     if (rn == 1)
  276       rp[1] = 0;
  277 #if GMP_NUMB_BITS > 32
  278     rp[1] = rp[1] << (GMP_NUMB_BITS - 32);
  279     rp[1] |= rp[0] >> 32;
  280     rp[0] &= 4294967295UL;
  281 #endif
  282     if (case_i)
  283       {  /* s < 2^53: case i) */
  284         x.s.exp = exp << 1;
  285         x.s.manl = rp[0];           /* 32 bits */
  286         x.s.manh = rp[1] & 1048575; /* 20 low bits */
  287         x.s.exp |= rp[1] >> 20;     /* 1 bit */
  288       }
  289     else /* s >= 2^53: case ii) */
  290       {
  291         x.s.exp = 1536 | (exp >> 1);
  292         x.s.manl = rp[0];
  293         x.s.manh = (rp[1] ^ 2097152) | ((exp & 1) << 19);
  294       }
  295   }
  296 #endif /* DPD_FORMAT */
  297   y.d = x.d;
  298   return y.d64;
  299 }
  300 #else
  301 /* portable version */
  302 static _Decimal64
  303 string_to_Decimal64 (char *s)
  304 {
  305   long int exp = 0;
  306   char m[17];
  307   long n = 0; /* mantissa length */
  308   char *endptr[1];
  309   _Decimal64 x = 0.0;
  310   int sign = 0;
  311 #ifdef DPD_FORMAT
  312   unsigned int G, d1, d2, d3, d4, d5;
  313 #endif
  314 
  315   /* read sign */
  316   if (*s == '-')
  317     {
  318       sign = 1;
  319       s ++;
  320     }
  321   /* read mantissa */
  322   while (ISDIGIT (*s))
  323     m[n++] = *s++;
  324   exp = n;
  325   if (*s == '.')
  326     {
  327       s ++;
  328       while (ISDIGIT (*s))
  329         m[n++] = *s++;
  330     }
  331   /* we have exp digits before decimal point, and a total of n digits */
  332   exp -= n; /* we will consider an integer mantissa */
  333   MPFR_ASSERTN(n <= 16);
  334   if (*s == 'E' || *s == 'e')
  335     exp += strtol (s + 1, endptr, 10);
  336   else
  337     *endptr = s;
  338   MPFR_ASSERTN(**endptr == '\0');
  339   MPFR_ASSERTN(-398 <= exp && exp <= (long) (385 - n));
  340   while (n < 16)
  341     {
  342       m[n++] = '0';
  343       exp --;
  344     }
  345   /* now n=16 and -398 <= exp <= 369 */
  346   m[n] = '\0';
  347 
  348   /* the number to convert is m[] * 10^exp where the mantissa is a 16-digit
  349      integer */
  350 
  351   /* compute biased exponent */
  352   exp += 398;
  353 
  354   MPFR_ASSERTN(exp >= -15);
  355   if (exp < 0)
  356     {
  357       int i;
  358       n = -exp;
  359       /* check the last n digits of the mantissa are zero */
  360       for (i = 1; i <= n; i++)
  361         MPFR_ASSERTN(m[16 - n] == '0');
  362       /* shift the first (16-n) digits to the right */
  363       for (i = 16 - n - 1; i >= 0; i--)
  364         m[i + n] = m[i];
  365       /* zero the first n digits */
  366       for (i = 0; i < n; i ++)
  367         m[i] = '0';
  368       exp = 0;
  369     }
  370 
  371   /* the number to convert is m[] * 10^(exp-398) */
  372   exp -= 398;
  373 
  374   for (n = 0; n < 16; n++)
  375     x = (_Decimal64) 10.0 * x + (_Decimal64) (m[n] - '0');
  376 
  377   /* multiply by 10^exp */
  378   if (exp > 0)
  379     {
  380       _Decimal64 ten16 = (double) 1e16; /* 10^16 is exactly representable
  381                                            in binary64 */
  382       _Decimal64 ten32 = ten16 * ten16;
  383       _Decimal64 ten64 = ten32 * ten32;
  384       _Decimal64 ten128 = ten64 * ten64;
  385       _Decimal64 ten256 = ten128 * ten128;
  386       if (exp >= 256)
  387         {
  388           x *= ten256;
  389           exp -= 256;
  390         }
  391       if (exp >= 128)
  392         {
  393           x *= ten128;
  394           exp -= 128;
  395         }
  396       if (exp >= 64)
  397         {
  398           x *= ten64;
  399           exp -= 64;
  400         }
  401       if (exp >= 32)
  402         {
  403           x *= ten32;
  404           exp -= 32;
  405         }
  406       if (exp >= 16)
  407         {
  408           x *= (_Decimal64) 10000000000000000.0;
  409           exp -= 16;
  410         }
  411       if (exp >= 8)
  412         {
  413           x *= (_Decimal64) 100000000.0;
  414           exp -= 8;
  415         }
  416       if (exp >= 4)
  417         {
  418           x *= (_Decimal64) 10000.0;
  419           exp -= 4;
  420         }
  421       if (exp >= 2)
  422         {
  423           x *= (_Decimal64) 100.0;
  424           exp -= 2;
  425         }
  426       if (exp >= 1)
  427         {
  428           x *= (_Decimal64) 10.0;
  429           exp -= 1;
  430         }
  431     }
  432   else if (exp < 0)
  433     {
  434       _Decimal64 ten16 = (double) 1e16; /* 10^16 is exactly representable
  435                                            in binary64 */
  436       _Decimal64 ten32 = ten16 * ten16;
  437       _Decimal64 ten64 = ten32 * ten32;
  438       _Decimal64 ten128 = ten64 * ten64;
  439       _Decimal64 ten256 = ten128 * ten128;
  440       if (exp <= -256)
  441         {
  442           x /= ten256;
  443           exp += 256;
  444         }
  445       if (exp <= -128)
  446         {
  447           x /= ten128;
  448           exp += 128;
  449         }
  450       if (exp <= -64)
  451         {
  452           x /= ten64;
  453           exp += 64;
  454         }
  455       if (exp <= -32)
  456         {
  457           x /= ten32;
  458           exp += 32;
  459         }
  460       if (exp <= -16)
  461         {
  462           x /= (_Decimal64) 10000000000000000.0;
  463           exp += 16;
  464         }
  465       if (exp <= -8)
  466         {
  467           x /= (_Decimal64) 100000000.0;
  468           exp += 8;
  469         }
  470       if (exp <= -4)
  471         {
  472           x /= (_Decimal64) 10000.0;
  473           exp += 4;
  474         }
  475       if (exp <= -2)
  476         {
  477           x /= (_Decimal64) 100.0;
  478           exp += 2;
  479         }
  480       if (exp <= -1)
  481         {
  482           x /= (_Decimal64) 10.0;
  483           exp += 1;
  484         }
  485     }
  486 
  487   if (sign)
  488     x = -x;
  489 
  490   return x;
  491 }
  492 #endif
  493 
  494 _Decimal64
  495 mpfr_get_decimal64 (mpfr_srcptr src, mpfr_rnd_t rnd_mode)
  496 {
  497   int negative;
  498   mpfr_exp_t e;
  499 
  500   /* the encoding of NaN, Inf, zero is the same under DPD or BID */
  501   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
  502     {
  503       if (MPFR_IS_NAN (src))
  504         return get_decimal64_nan ();
  505 
  506       negative = MPFR_IS_NEG (src);
  507 
  508       if (MPFR_IS_INF (src))
  509         return get_decimal64_inf (negative);
  510 
  511       MPFR_ASSERTD (MPFR_IS_ZERO(src));
  512       return get_decimal64_zero (negative);
  513     }
  514 
  515   e = MPFR_GET_EXP (src);
  516   negative = MPFR_IS_NEG (src);
  517 
  518   if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDA))
  519     rnd_mode = negative ? MPFR_RNDD : MPFR_RNDU;
  520 
  521   /* the smallest decimal64 number is 10^(-398),
  522      with 2^(-1323) < 10^(-398) < 2^(-1322) */
  523   if (MPFR_UNLIKELY (e < -1323)) /* src <= 2^(-1324) < 1/2*10^(-398) */
  524     {
  525       if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDN
  526           || (rnd_mode == MPFR_RNDD && negative == 0)
  527           || (rnd_mode == MPFR_RNDU && negative != 0))
  528         return get_decimal64_zero (negative);
  529       else /* return the smallest non-zero number */
  530         return get_decimal64_min (negative);
  531     }
  532   /* the largest decimal64 number is just below 10^(385) < 2^1279 */
  533   else if (MPFR_UNLIKELY (e > 1279)) /* then src >= 2^1279 */
  534     {
  535       if (rnd_mode == MPFR_RNDZ
  536           || (rnd_mode == MPFR_RNDU && negative != 0)
  537           || (rnd_mode == MPFR_RNDD && negative == 0))
  538         return get_decimal64_max (negative);
  539       else
  540         return get_decimal64_inf (negative);
  541     }
  542   else
  543     {
  544       /* we need to store the sign (1), the mantissa (16), and the terminating
  545          character, thus we need at least 18 characters in s */
  546       char s[23];
  547       mpfr_get_str (s, &e, 10, 16, src, rnd_mode);
  548       /* the smallest normal number is 1.000...000E-383,
  549          which corresponds to s=[0.]1000...000 and e=-382 */
  550       if (e < -382)
  551         {
  552           /* the smallest subnormal number is 0.000...001E-383 = 1E-398,
  553              which corresponds to s=[0.]1000...000 and e=-397 */
  554           if (e < -397)
  555             {
  556               if (rnd_mode == MPFR_RNDN && e == -398)
  557                 {
  558                   /* If 0.5E-398 < |src| < 1E-398 (smallest subnormal),
  559                      src should round to +/- 1E-398 in MPFR_RNDN. */
  560                   mpfr_get_str (s, &e, 10, 1, src, MPFR_RNDA);
  561                   return e == -398 && s[negative] <= '5' ?
  562                     get_decimal64_zero (negative) :
  563                     get_decimal64_min (negative);
  564                 }
  565               if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDN
  566                   || (rnd_mode == MPFR_RNDD && negative == 0)
  567                   || (rnd_mode == MPFR_RNDU && negative != 0))
  568                 return get_decimal64_zero (negative);
  569               else /* return the smallest non-zero number */
  570                 return get_decimal64_min (negative);
  571             }
  572           else
  573             {
  574               mpfr_exp_t e2;
  575               long digits = 16 - (-382 - e);
  576               /* if e = -397 then 16 - (-382 - e) = 1 */
  577               mpfr_get_str (s, &e2, 10, digits, src, rnd_mode);
  578               /* Warning: we can have e2 = e + 1 here, when rounding to
  579                  nearest or away from zero. */
  580               s[negative + digits] = 'E';
  581               sprintf (s + negative + digits + 1, "%ld",
  582                        (long int)e2 - digits);
  583               return string_to_Decimal64 (s);
  584             }
  585         }
  586       /* the largest number is 9.999...999E+384,
  587          which corresponds to s=[0.]9999...999 and e=385 */
  588       else if (e > 385)
  589         {
  590           if (rnd_mode == MPFR_RNDZ
  591               || (rnd_mode == MPFR_RNDU && negative != 0)
  592               || (rnd_mode == MPFR_RNDD && negative == 0))
  593             return get_decimal64_max (negative);
  594           else
  595             return get_decimal64_inf (negative);
  596         }
  597       else /* -382 <= e <= 385 */
  598         {
  599           s[16 + negative] = 'E';
  600           sprintf (s + 17 + negative, "%ld", (long int)e - 16);
  601           return string_to_Decimal64 (s);
  602         }
  603     }
  604 }
  605 
  606 #endif /* MPFR_WANT_DECIMAL_FLOATS */