xref: /netbsd-src/external/lgpl3/mpfr/dist/src/get_d64.c (revision e670fd5c413e99c2f6a37901bb21c537fcd322d2)
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 <http://www.open-std.org/jtc1/sc22/wg14/www/projects#24732>.
7 
8 Copyright 2006-2020 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.0dd : 0.0dd;
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         return get_decimal64_nan ();
517 
518       negative = MPFR_IS_NEG (src);
519 
520       if (MPFR_IS_INF (src))
521         return get_decimal64_inf (negative);
522 
523       MPFR_ASSERTD (MPFR_IS_ZERO(src));
524       return get_decimal64_zero (negative);
525     }
526 
527   e = MPFR_GET_EXP (src);
528   negative = MPFR_IS_NEG (src);
529 
530   MPFR_UPDATE2_RND_MODE (rnd_mode, MPFR_SIGN (src));
531 
532   /* now rnd_mode is RNDN, RNDF, RNDA or RNDZ */
533 
534   /* the smallest decimal64 number is 10^(-398),
535      with 2^(-1323) < 10^(-398) < 2^(-1322) */
536   if (MPFR_UNLIKELY (e < -1323)) /* src <= 2^(-1324) < 1/2*10^(-398) */
537     {
538       if (rnd_mode != MPFR_RNDA)
539         return get_decimal64_zero (negative);
540       else /* RNDA: return the smallest non-zero number */
541         return get_decimal64_min (negative);
542     }
543   /* the largest decimal64 number is just below 10^385 < 2^1279 */
544   else if (MPFR_UNLIKELY (e > 1279)) /* then src >= 2^1279 */
545     {
546       if (rnd_mode == MPFR_RNDZ)
547         return get_decimal64_max (negative);
548       else /* RNDN, RNDA, RNDF: round away */
549         return get_decimal64_inf (negative);
550     }
551   else
552     {
553       /* we need to store the sign (1 character), the significand (at most 16
554          characters), the exponent part (at most 5 characters for "E-398"),
555          and the terminating character, thus we need at least 23 characters */
556       char s[23];
557       mpfr_get_str (s, &e, 10, 16, src, rnd_mode);
558       /* the smallest normal number is 1.000...000E-383,
559          which corresponds to s=[0.]1000...000 and e=-382 */
560       if (e < -382)
561         {
562           /* the smallest subnormal number is 0.000...001E-383 = 1E-398,
563              which corresponds to s=[0.]1000...000 and e=-397 */
564           if (e < -397)
565             {
566               if (rnd_mode == MPFR_RNDN && e == -398)
567                 {
568                   /* If 0.5E-398 < |src| < 1E-398 (smallest subnormal),
569                      src should round to +/- 1E-398 in MPFR_RNDN. */
570                   mpfr_get_str (s, &e, 10, 1, src, MPFR_RNDA);
571                   return e == -398 && s[negative] <= '5' ?
572                     get_decimal64_zero (negative) :
573                     get_decimal64_min (negative);
574                 }
575               if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDN)
576                 return get_decimal64_zero (negative);
577               else /* RNDA or RNDF: return the smallest non-zero number */
578                 return get_decimal64_min (negative);
579             }
580           else
581             {
582               mpfr_exp_t e2;
583               long digits = 16 - (-382 - e);
584               /* if e = -397 then 16 - (-382 - e) = 1 */
585               mpfr_get_str (s, &e2, 10, digits, src, rnd_mode);
586               /* Warning: we can have e2 = e + 1 here, when rounding to
587                  nearest or away from zero. */
588               s[negative + digits] = 'E';
589               sprintf (s + negative + digits + 1, "%ld",
590                        (long int)e2 - digits);
591               return string_to_Decimal64 (s);
592             }
593         }
594       /* the largest number is 9.999...999E+384,
595          which corresponds to s=[0.]9999...999 and e=385 */
596       else if (e > 385)
597         {
598           if (rnd_mode == MPFR_RNDZ)
599             return get_decimal64_max (negative);
600           else /* RNDN, RNDA, RNDF: round away */
601             return get_decimal64_inf (negative);
602         }
603       else /* -382 <= e <= 385 */
604         {
605           s[16 + negative] = 'E';
606           sprintf (s + 17 + negative, "%ld", (long int)e - 16);
607           return string_to_Decimal64 (s);
608         }
609     }
610 }
611 
612 #endif /* MPFR_WANT_DECIMAL_FLOATS */
613