1/* mpfr_set_decimal128 -- convert a IEEE 754r decimal128 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 Caramel 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#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 */
39static 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*/
126static void
127decimal128_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
257static void
258decimal128_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 */
481int
482mpfr_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 */
497