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