1/* mpfr_set_decimal64 -- convert an IEEE 754-2008 decimal64 float to 2 a multiple precision floating-point number 3 4See https://gcc.gnu.org/legacy-ml/gcc/2006-06/msg00691.html, 5https://gcc.gnu.org/onlinedocs/gcc/Decimal-Float.html, 6and TR 24732 <https://www.open-std.org/jtc1/sc22/wg14/www/projects#24732>. 7 8Copyright 2006-2023 Free Software Foundation, Inc. 9Contributed by the AriC and Caramba 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 25https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2651 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 */ 35static 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 */ 115static void 116decimal64_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 284static void 285decimal64_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 /* we don't propagate the sign bit */ 293 sprintf (s, "NaN"); /* sprintf puts a final '\0' */ 294 return; 295 } 296 else if (MPFR_UNLIKELY (d > DEC64_MAX)) /* +Inf */ 297 { 298 sprintf (s, "Inf"); 299 return; 300 } 301 else if (MPFR_UNLIKELY (d < -DEC64_MAX)) /* -Inf */ 302 { 303 sprintf (s, "-Inf"); 304 return; 305 } 306 307 /* now d is neither NaN nor +Inf nor -Inf */ 308 309 if (d < (_Decimal64) 0.0) 310 { 311 sign = 1; 312 d = -d; 313 } 314 else if (d == (_Decimal64) -0.0) 315 { /* Warning: the comparison d == -0.0 returns true for d = 0.0 too, 316 copy code from set_d.c here. We first compare to the +0.0 bitstring, 317 in case +0.0 and -0.0 are represented identically. */ 318 double dd = (double) d, poszero = +0.0, negzero = DBL_NEG_ZERO; 319 if (memcmp (&dd, &poszero, sizeof(double)) != 0 && 320 memcmp (&dd, &negzero, sizeof(double)) == 0) 321 { 322 sign = 1; 323 d = -d; 324 } 325 } 326 327 /* now normalize d in [0.1, 1[ */ 328 if (d >= (_Decimal64) 1.0) 329 { 330 _Decimal64 ten16 = (double) 1e16; /* 10^16 is exactly representable 331 in binary64 */ 332 _Decimal64 ten32 = ten16 * ten16; 333 _Decimal64 ten64 = ten32 * ten32; 334 _Decimal64 ten128 = ten64 * ten64; 335 _Decimal64 ten256 = ten128 * ten128; 336 337 if (d >= ten256) 338 { 339 d /= ten256; 340 exp += 256; 341 } 342 if (d >= ten128) 343 { 344 d /= ten128; 345 exp += 128; 346 } 347 if (d >= ten64) 348 { 349 d /= ten64; 350 exp += 64; 351 } 352 if (d >= ten32) 353 { 354 d /= ten32; 355 exp += 32; 356 } 357 if (d >= (_Decimal64) 10000000000000000.0) 358 { 359 d /= (_Decimal64) 10000000000000000.0; 360 exp += 16; 361 } 362 if (d >= (_Decimal64) 100000000.0) 363 { 364 d /= (_Decimal64) 100000000.0; 365 exp += 8; 366 } 367 if (d >= (_Decimal64) 10000.0) 368 { 369 d /= (_Decimal64) 10000.0; 370 exp += 4; 371 } 372 if (d >= (_Decimal64) 100.0) 373 { 374 d /= (_Decimal64) 100.0; 375 exp += 2; 376 } 377 if (d >= (_Decimal64) 10.0) 378 { 379 d /= (_Decimal64) 10.0; 380 exp += 1; 381 } 382 if (d >= (_Decimal64) 1.0) 383 { 384 d /= (_Decimal64) 10.0; 385 exp += 1; 386 } 387 } 388 else /* d < 1.0 */ 389 { 390 _Decimal64 ten16, ten32, ten64, ten128, ten256; 391 392 ten16 = (double) 1e16; /* 10^16 is exactly representable in binary64 */ 393 ten16 = (_Decimal64) 1.0 / ten16; /* 10^(-16), exact */ 394 ten32 = ten16 * ten16; 395 ten64 = ten32 * ten32; 396 ten128 = ten64 * ten64; 397 ten256 = ten128 * ten128; 398 399 if (d < ten256) 400 { 401 d /= ten256; 402 exp -= 256; 403 } 404 if (d < ten128) 405 { 406 d /= ten128; 407 exp -= 128; 408 } 409 if (d < ten64) 410 { 411 d /= ten64; 412 exp -= 64; 413 } 414 if (d < ten32) 415 { 416 d /= ten32; 417 exp -= 32; 418 } 419 /* the double constant 0.0000000000000001 is 2028240960365167/2^104, 420 which should be rounded to 1e-16 in _Decimal64 */ 421 if (d < (_Decimal64) 0.0000000000000001) 422 { 423 d *= (_Decimal64) 10000000000000000.0; 424 exp -= 16; 425 } 426 /* the double constant 0.00000001 is 3022314549036573/2^78, 427 which should be rounded to 1e-8 in _Decimal64 */ 428 if (d < (_Decimal64) 0.00000001) 429 { 430 d *= (_Decimal64) 100000000.0; 431 exp -= 8; 432 } 433 /* the double constant 0.0001 is 7378697629483821/2^66, 434 which should be rounded to 1e-4 in _Decimal64 */ 435 if (d < (_Decimal64) 0.0001) 436 { 437 d *= (_Decimal64) 10000.0; 438 exp -= 4; 439 } 440 /* the double constant 0.01 is 5764607523034235/2^59, 441 which should be rounded to 1e-2 in _Decimal64 */ 442 if (d < (_Decimal64) 0.01) 443 { 444 d *= (_Decimal64) 100.0; 445 exp -= 2; 446 } 447 /* the double constant 0.1 is 3602879701896397/2^55, 448 which should be rounded to 1e-1 in _Decimal64 */ 449 if (d < (_Decimal64) 0.1) 450 { 451 d *= (_Decimal64) 10.0; 452 exp -= 1; 453 } 454 } 455 456 /* now 0.1 <= d < 1 */ 457 if (sign == 1) 458 *s++ = '-'; 459 *s++ = '0'; 460 *s++ = '.'; 461 for (n = 0; n < 16; n++) 462 { 463 double e; 464 int r; 465 466 d *= (_Decimal64) 10.0; 467 e = (double) d; 468 r = (int) e; 469 *s++ = '0' + r; 470 d -= (_Decimal64) r; 471 } 472 MPFR_ASSERTN(d == (_Decimal64) 0.0); 473 if (exp != 0) 474 sprintf (s, "E%d", exp); /* adds a final '\0' */ 475 else 476 *s = '\0'; 477} 478 479#endif /* definition of decimal64_to_string (DPD, BID, or portable) */ 480 481/* the IEEE754-2008 decimal64 format has 16 digits, with emax=384, 482 emin=1-emax=-383 */ 483int 484mpfr_set_decimal64 (mpfr_ptr r, _Decimal64 d, mpfr_rnd_t rnd_mode) 485{ 486 char s[25]; /* need 1 character for sign, 487 2 characters for '0.' 488 16 characters for significand, 489 1 character for exponent 'E', 490 4 characters for exponent (including sign), 491 1 character for terminating \0. */ 492 493 decimal64_to_string (s, d); 494 MPFR_LOG_MSG (("string: %s\n", s)); 495 return mpfr_strtofr (r, s, NULL, 10, rnd_mode); 496} 497 498#endif /* MPFR_WANT_DECIMAL_FLOATS */ 499