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