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 <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 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
decimal128_to_string(char * s,_Decimal128 d)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
decimal128_to_string(char * s,_Decimal128 d)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 sprintf (s, "NaN"); /* sprintf puts a final '\0' */
282 return;
283 }
284 else if (MPFR_UNLIKELY (d > DEC128_MAX)) /* +Inf */
285 {
286 sprintf (s, "Inf");
287 return;
288 }
289 else if (MPFR_UNLIKELY (d < -DEC128_MAX)) /* -Inf */
290 {
291 sprintf (s, "-Inf");
292 return;
293 }
294
295 /* now d is neither NaN nor +Inf nor -Inf */
296
297 if (d < (_Decimal128) 0.0)
298 {
299 sign = 1;
300 d = -d;
301 }
302 else if (d == (_Decimal128) -0.0)
303 { /* Warning: the comparison d == -0.0 returns true for d = 0.0 too,
304 copy code from set_d.c here. We first compare to the +0.0 bitstring,
305 in case +0.0 and -0.0 are represented identically. */
306 double dd = (double) d, poszero = +0.0, negzero = DBL_NEG_ZERO;
307 if (memcmp (&dd, &poszero, sizeof(double)) != 0 &&
308 memcmp (&dd, &negzero, sizeof(double)) == 0)
309 {
310 sign = 1;
311 d = -d;
312 }
313 }
314
315 /* now normalize d in [0.1, 1[ */
316 if (d >= (_Decimal128) 1.0)
317 {
318 if (d >= ten4096)
319 {
320 d /= ten4096;
321 exp += 4096;
322 }
323 if (d >= ten2048)
324 {
325 d /= ten2048;
326 exp += 2048;
327 }
328 if (d >= ten1024)
329 {
330 d /= ten1024;
331 exp += 1024;
332 }
333 if (d >= ten512)
334 {
335 d /= ten512;
336 exp += 512;
337 }
338 if (d >= ten256)
339 {
340 d /= ten256;
341 exp += 256;
342 }
343 if (d >= ten128)
344 {
345 d /= ten128;
346 exp += 128;
347 }
348 if (d >= ten64)
349 {
350 d /= ten64;
351 exp += 64;
352 }
353 if (d >= ten32)
354 {
355 d /= ten32;
356 exp += 32;
357 }
358 if (d >= ten16)
359 {
360 d /= ten16;
361 exp += 16;
362 }
363 if (d >= ten8)
364 {
365 d /= ten8;
366 exp += 8;
367 }
368 if (d >= ten4)
369 {
370 d /= ten4;
371 exp += 4;
372 }
373 if (d >= ten2)
374 {
375 d /= ten2;
376 exp += 2;
377 }
378 while (d >= 1)
379 {
380 d /= ten;
381 exp += 1;
382 }
383 }
384 else /* d < 1.0 */
385 {
386 if (d < 1 / ten4096)
387 {
388 d *= ten4096;
389 exp -= 4096;
390 }
391 if (d < 1 / ten2048)
392 {
393 d *= ten2048;
394 exp -= 2048;
395 }
396 if (d < 1 / ten1024)
397 {
398 d *= ten1024;
399 exp -= 1024;
400 }
401 if (d < 1 / ten512)
402 {
403 d *= ten512;
404 exp -= 512;
405 }
406 if (d < 1 / ten256)
407 {
408 d *= ten256;
409 exp -= 256;
410 }
411 if (d < 1 / ten128)
412 {
413 d *= ten128;
414 exp -= 128;
415 }
416 if (d < 1 / ten64)
417 {
418 d *= ten64;
419 exp -= 64;
420 }
421 if (d < 1 / ten32)
422 {
423 d *= ten32;
424 exp -= 32;
425 }
426 if (d < 1 / ten16)
427 {
428 d *= ten16;
429 exp -= 16;
430 }
431 if (d < 1 / ten8)
432 {
433 d *= ten8;
434 exp -= 8;
435 }
436 if (d < 1 / ten4)
437 {
438 d *= ten4;
439 exp -= 4;
440 }
441 if (d < 1 / ten2)
442 {
443 d *= ten2;
444 exp -= 2;
445 }
446 if (d < 1 / ten)
447 {
448 d *= ten;
449 exp -= 1;
450 }
451 }
452
453 /* now 0.1 <= d < 1 */
454 if (sign == 1)
455 *s++ = '-';
456 *s++ = '0';
457 *s++ = '.';
458 for (n = 0; n < 34; n++)
459 {
460 int r;
461
462 MPFR_ASSERTD (d < 1);
463 d *= 10;
464 MPFR_ASSERTD (d < 10);
465 r = (int) d;
466 *s++ = '0' + r;
467 d -= (_Decimal128) r;
468 }
469 MPFR_ASSERTN (d == 0);
470 if (exp != 0)
471 sprintf (s, "E%d", exp); /* adds a final '\0' */
472 else
473 *s = '\0';
474 }
475
476 #endif /* definition of decimal128_to_string (DPD, BID, or portable) */
477
478 /* the IEEE754-2008 decimal128 format has 34 digits, with emax=6144,
479 emin=1-emax=-6143 */
480 int
mpfr_set_decimal128(mpfr_ptr r,_Decimal128 d,mpfr_rnd_t rnd_mode)481 mpfr_set_decimal128 (mpfr_ptr r, _Decimal128 d, mpfr_rnd_t rnd_mode)
482 {
483 char s[44]; /* need 1 character for sign,
484 2 characters for '0.'
485 34 characters for significand,
486 1 character for exponent 'E',
487 5 characters for exponent (including sign),
488 1 character for terminating \0. */
489
490 decimal128_to_string (s, d);
491 MPFR_LOG_MSG (("string: %s\n", s));
492 return mpfr_strtofr (r, s, NULL, 10, rnd_mode);
493 }
494
495 #endif /* MPFR_WANT_DECIMAL_FLOATS */
496