1/* mpfr_get_decimal64 -- convert a multiple precision floating-point number 2 to a IEEE 754r decimal64 float 3 4See http://gcc.gnu.org/ml/gcc/2006-06/msg00691.html, 5http://gcc.gnu.org/onlinedocs/gcc/Decimal-Float.html, 6and TR 24732 <http://www.open-std.org/jtc1/sc22/wg14/www/projects#24732>. 7 8Copyright 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 9Contributed by the Arenaire and Cacao projects, INRIA. 10 11This file is part of the GNU MPFR Library. 12 13The GNU MPFR Library is free software; you can redistribute it and/or modify 14it under the terms of the GNU Lesser General Public License as published by 15the Free Software Foundation; either version 3 of the License, or (at your 16option) any later version. 17 18The GNU MPFR Library is distributed in the hope that it will be useful, but 19WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 20or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 21License for more details. 22 23You should have received a copy of the GNU Lesser General Public License 24along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 25http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2651 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#ifdef DPD_FORMAT 36static 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}; 103#endif 104 105/* construct a decimal64 NaN */ 106static _Decimal64 107get_decimal64_nan (void) 108{ 109 union ieee_double_extract x; 110 union ieee_double_decimal64 y; 111 112 x.s.exp = 1984; /* G[0]..G[4] = 11111: quiet NaN */ 113 y.d = x.d; 114 return y.d64; 115} 116 117/* construct the decimal64 Inf with given sign */ 118static _Decimal64 119get_decimal64_inf (int negative) 120{ 121 union ieee_double_extract x; 122 union ieee_double_decimal64 y; 123 124 x.s.sig = (negative) ? 1 : 0; 125 x.s.exp = 1920; /* G[0]..G[4] = 11110: Inf */ 126 y.d = x.d; 127 return y.d64; 128} 129 130/* construct the decimal64 zero with given sign */ 131static _Decimal64 132get_decimal64_zero (int negative) 133{ 134 union ieee_double_decimal64 y; 135 136 /* zero has the same representation in binary64 and decimal64 */ 137 y.d = negative ? DBL_NEG_ZERO : 0.0; 138 return y.d64; 139} 140 141/* construct the decimal64 smallest non-zero with given sign */ 142static _Decimal64 143get_decimal64_min (int negative) 144{ 145 union ieee_double_extract x; 146 147 x.s.sig = (negative) ? 1 : 0; 148 x.s.exp = 0; 149 x.s.manh = 0; 150 x.s.manl = 1; 151 return x.d; 152} 153 154/* construct the decimal64 largest finite number with given sign */ 155static _Decimal64 156get_decimal64_max (int negative) 157{ 158 union ieee_double_extract x; 159 160 x.s.sig = (negative) ? 1 : 0; 161 x.s.exp = 1919; 162 x.s.manh = 1048575; /* 2^20-1 */ 163 x.s.manl = ~0; 164 return x.d; 165} 166 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. 174*/ 175static _Decimal64 176string_to_Decimal64 (char *s) 177{ 178 long int exp = 0; 179 char m[17]; 180 long n = 0; /* mantissa length */ 181 char *endptr[1]; 182 union ieee_double_extract x; 183 union ieee_double_decimal64 y; 184#ifdef DPD_FORMAT 185 unsigned int G, d1, d2, d3, d4, d5; 186#endif 187 188 /* read sign */ 189 if (*s == '-') 190 { 191 x.s.sig = 1; 192 s ++; 193 } 194 else 195 x.s.sig = 0; 196 /* read mantissa */ 197 while (ISDIGIT (*s)) 198 m[n++] = *s++; 199 exp = n; 200 if (*s == '.') 201 { 202 s ++; 203 while (ISDIGIT (*s)) 204 m[n++] = *s++; 205 } 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); 211 else 212 *endptr = s; 213 MPFR_ASSERTN(**endptr == '\0'); 214 MPFR_ASSERTN(-398 <= exp && exp <= (long) (385 - n)); 215 while (n < 16) 216 { 217 m[n++] = '0'; 218 exp --; 219 } 220 /* now n=16 and -398 <= exp <= 369 */ 221 m[n] = '\0'; 222 223 /* compute biased exponent */ 224 exp += 398; 225 226 MPFR_ASSERTN(exp >= -15); 227 if (exp < 0) 228 { 229 int i; 230 n = -exp; 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--) 236 m[i + n] = m[i]; 237 /* zero the first n digits */ 238 for (i = 0; i < n; i ++) 239 m[i] = '0'; 240 exp = 0; 241 } 242 243 /* now convert to DPD or BID */ 244#ifdef DPD_FORMAT 245#define CH(d) (d - '0') 246 if (m[0] >= '8') 247 G = (3 << 11) | ((exp & 768) << 1) | ((CH(m[0]) & 1) << 8); 248 else 249 G = ((exp & 768) << 3) | (CH(m[0]) << 8); 250 /* now the most 5 significant bits of G are filled */ 251 G |= exp & 255; 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 */ 257 x.s.exp = G >> 2; 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 */ 262 { 263 mp_size_t rn; 264 mp_limb_t rp[2]; 265 int case_i = strcmp (m, "9007199254740992") < 0; 266 267 for (n = 0; n < 16; n++) 268 m[n] -= '0'; 269 rn = mpn_set_str (rp, (unsigned char *) m, 16, 10); 270 if (rn == 1) 271 rp[1] = 0; 272#if GMP_NUMB_BITS > 32 273 rp[1] = rp[1] << (GMP_NUMB_BITS - 32); 274 rp[1] |= rp[0] >> 32; 275 rp[0] &= 4294967295UL; 276#endif 277 if (case_i) 278 { /* s < 2^53: case i) */ 279 x.s.exp = exp << 1; 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 */ 283 } 284 else /* s >= 2^53: case ii) */ 285 { 286 x.s.exp = 1536 | (exp >> 1); 287 x.s.manl = rp[0]; 288 x.s.manh = (rp[1] ^ 2097152) | ((exp & 1) << 19); 289 } 290 } 291#endif /* DPD_FORMAT */ 292 y.d = x.d; 293 return y.d64; 294} 295 296_Decimal64 297mpfr_get_decimal64 (mpfr_srcptr src, mpfr_rnd_t rnd_mode) 298{ 299 int negative; 300 mpfr_exp_t e; 301 302 /* the encoding of NaN, Inf, zero is the same under DPD or BID */ 303 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src))) 304 { 305 if (MPFR_IS_NAN (src)) 306 return get_decimal64_nan (); 307 308 negative = MPFR_IS_NEG (src); 309 310 if (MPFR_IS_INF (src)) 311 return get_decimal64_inf (negative); 312 313 MPFR_ASSERTD (MPFR_IS_ZERO(src)); 314 return get_decimal64_zero (negative); 315 } 316 317 e = MPFR_GET_EXP (src); 318 negative = MPFR_IS_NEG (src); 319 320 if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDA)) 321 rnd_mode = negative ? MPFR_RNDD : MPFR_RNDU; 322 323 /* the smallest decimal64 number is 10^(-398), 324 with 2^(-1323) < 10^(-398) < 2^(-1322) */ 325 if (MPFR_UNLIKELY (e < -1323)) /* src <= 2^(-1324) < 1/2*10^(-398) */ 326 { 327 if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDN 328 || (rnd_mode == MPFR_RNDD && negative == 0) 329 || (rnd_mode == MPFR_RNDU && negative != 0)) 330 return get_decimal64_zero (negative); 331 else /* return the smallest non-zero number */ 332 return get_decimal64_min (negative); 333 } 334 /* the largest decimal64 number is just below 10^(385) < 2^1279 */ 335 else if (MPFR_UNLIKELY (e > 1279)) /* then src >= 2^1279 */ 336 { 337 if (MPFR_RNDZ || (rnd_mode == MPFR_RNDU && negative != 0) 338 || (rnd_mode == MPFR_RNDD && negative == 0)) 339 return get_decimal64_max (negative); 340 else 341 return get_decimal64_inf (negative); 342 } 343 else 344 { 345 /* we need to store the sign (1), the mantissa (16), and the terminating 346 character, thus we need at least 18 characters in s */ 347 char s[23]; 348 mpfr_get_str (s, &e, 10, 16, src, rnd_mode); 349 /* the smallest normal number is 1.000...000E-383, 350 which corresponds to s=[0.]1000...000 and e=-382 */ 351 if (e < -382) 352 { 353 /* the smallest subnormal number is 0.000...001E-383 = 1E-398, 354 which corresponds to s=[0.]1000...000 and e=-397 */ 355 if (e < -397) 356 { 357 if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDN 358 || (rnd_mode == MPFR_RNDD && negative == 0) 359 || (rnd_mode == MPFR_RNDU && negative != 0)) 360 return get_decimal64_zero (negative); 361 else /* return the smallest non-zero number */ 362 return get_decimal64_min (negative); 363 } 364 else 365 { 366 mpfr_exp_t e2; 367 long digits = 16 - (-382 - e); 368 /* if e = -397 then 16 - (-382 - e) = 1 */ 369 mpfr_get_str (s, &e2, 10, digits, src, rnd_mode); 370 /* Warning: we can have e2 = e + 1 here, when rounding to 371 nearest or away from zero. */ 372 s[negative + digits] = 'E'; 373 sprintf (s + negative + digits + 1, "%ld", 374 (long int)e2 - digits); 375 return string_to_Decimal64 (s); 376 } 377 } 378 /* the largest number is 9.999...999E+384, 379 which corresponds to s=[0.]9999...999 and e=385 */ 380 else if (e > 385) 381 { 382 if (MPFR_RNDZ || (rnd_mode == MPFR_RNDU && negative != 0) 383 || (rnd_mode == MPFR_RNDD && negative == 0)) 384 return get_decimal64_max (negative); 385 else 386 return get_decimal64_inf (negative); 387 } 388 else /* -382 <= e <= 385 */ 389 { 390 s[16 + negative] = 'E'; 391 sprintf (s + 17 + negative, "%ld", (long int)e - 16); 392 return string_to_Decimal64 (s); 393 } 394 } 395} 396 397#endif /* MPFR_WANT_DECIMAL_FLOATS */ 398