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