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