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