xref: /netbsd-src/external/lgpl3/mpfr/dist/src/get_d64.c (revision ba125506a622fe649968631a56eba5d42ff57863)
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
get_decimal64_nan(void)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
get_decimal64_inf(int negative)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
get_decimal64_zero(int negative)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
get_decimal64_min(int negative)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
get_decimal64_max(int negative)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
string_to_Decimal64(char * s)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
string_to_Decimal64(char * s)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
mpfr_get_decimal64(mpfr_srcptr src,mpfr_rnd_t rnd_mode)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 */
616