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 <https://www.open-std.org/jtc1/sc22/wg14/www/projects#24732>.
7
8 Copyright 2006-2023 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 /* 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 */
483 int
484 mpfr_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 */