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