source/libs/mpfr/mpfr-src/src/get_d64.c

   1 /* mpfr_get_decimal64 -- convert a multiple precision floating-point number
   2                          to a IEEE 754r decimal64 float
   3
   4 See http://gcc.gnu.org/ml/gcc/2006-06/msg00691.html,
   5 http://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-2016 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 http://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 #include <stdlib.h> /* for strtol */
  29 #include "mpfr-impl.h"
  30
  31 #define ISDIGIT(c) ('0' <= c && c <= '9')
  32
  33 #ifdef MPFR_WANT_DECIMAL_FLOATS
  34
  35 #ifndef DEC64_MAX
  36 # define DEC64_MAX 9.999999999999999E384dd
  37 #endif
  38
  39 #ifdef DPD_FORMAT
  40 static int T[1000] = {
  41   0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 32,
  42   33, 34, 35, 36, 37, 38, 39, 40, 41, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57,
  43   64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 80, 81, 82, 83, 84, 85, 86, 87, 88,
  44   89, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 112, 113, 114, 115, 116,
  45   117, 118, 119, 120, 121, 10, 11, 42, 43, 74, 75, 106, 107, 78, 79, 26, 27,
  46   58, 59, 90, 91, 122, 123, 94, 95, 128, 129, 130, 131, 132, 133, 134, 135,
  47   136, 137, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 160, 161, 162,
  48   163, 164, 165, 166, 167, 168, 169, 176, 177, 178, 179, 180, 181, 182, 183,
  49   184, 185, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 208, 209, 210,
  50   211, 212, 213, 214, 215, 216, 217, 224, 225, 226, 227, 228, 229, 230, 231,
  51   232, 233, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 138, 139, 170,
  52   171, 202, 203, 234, 235, 206, 207, 154, 155, 186, 187, 218, 219, 250, 251,
  53   222, 223, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 272, 273, 274,
  54   275, 276, 277, 278, 279, 280, 281, 288, 289, 290, 291, 292, 293, 294, 295,
  55   296, 297, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 320, 321, 322,
  56   323, 324, 325, 326, 327, 328, 329, 336, 337, 338, 339, 340, 341, 342, 343,
  57   344, 345, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 368, 369, 370,
  58   371, 372, 373, 374, 375, 376, 377, 266, 267, 298, 299, 330, 331, 362, 363,
  59   334, 335, 282, 283, 314, 315, 346, 347, 378, 379, 350, 351, 384, 385, 386,
  60   387, 388, 389, 390, 391, 392, 393, 400, 401, 402, 403, 404, 405, 406, 407,
  61   408, 409, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 432, 433, 434,
  62   435, 436, 437, 438, 439, 440, 441, 448, 449, 450, 451, 452, 453, 454, 455,
  63   456, 457, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 480, 481, 482,
  64   483, 484, 485, 486, 487, 488, 489, 496, 497, 498, 499, 500, 501, 502, 503,
  65   504, 505, 394, 395, 426, 427, 458, 459, 490, 491, 462, 463, 410, 411, 442,
  66   443, 474, 475, 506, 507, 478, 479, 512, 513, 514, 515, 516, 517, 518, 519,
  67   520, 521, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 544, 545, 546,
  68   547, 548, 549, 550, 551, 552, 553, 560, 561, 562, 563, 564, 565, 566, 567,
  69   568, 569, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 592, 593, 594,
  70   595, 596, 597, 598, 599, 600, 601, 608, 609, 610, 611, 612, 613, 614, 615,
  71   616, 617, 624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 522, 523, 554,
  72   555, 586, 587, 618, 619, 590, 591, 538, 539, 570, 571, 602, 603, 634, 635,
  73   606, 607, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 656, 657, 658,
  74   659, 660, 661, 662, 663, 664, 665, 672, 673, 674, 675, 676, 677, 678, 679,
  75   680, 681, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 704, 705, 706,
  76   707, 708, 709, 710, 711, 712, 713, 720, 721, 722, 723, 724, 725, 726, 727,
  77   728, 729, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 752, 753, 754,
  78   755, 756, 757, 758, 759, 760, 761, 650, 651, 682, 683, 714, 715, 746, 747,
  79   718, 719, 666, 667, 698, 699, 730, 731, 762, 763, 734, 735, 768, 769, 770,
  80   771, 772, 773, 774, 775, 776, 777, 784, 785, 786, 787, 788, 789, 790, 791,
  81   792, 793, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 816, 817, 818,
  82   819, 820, 821, 822, 823, 824, 825, 832, 833, 834, 835, 836, 837, 838, 839,
  83   840, 841, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 864, 865, 866,
  84   867, 868, 869, 870, 871, 872, 873, 880, 881, 882, 883, 884, 885, 886, 887,
  85   888, 889, 778, 779, 810, 811, 842, 843, 874, 875, 846, 847, 794, 795, 826,
  86   827, 858, 859, 890, 891, 862, 863, 896, 897, 898, 899, 900, 901, 902, 903,
  87   904, 905, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 928, 929, 930,
  88   931, 932, 933, 934, 935, 936, 937, 944, 945, 946, 947, 948, 949, 950, 951,
  89   952, 953, 960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 976, 977, 978,
  90   979, 980, 981, 982, 983, 984, 985, 992, 993, 994, 995, 996, 997, 998, 999,
  91   1000, 1001, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 906,
  92   907, 938, 939, 970, 971, 1002, 1003, 974, 975, 922, 923, 954, 955, 986,
  93   987, 1018, 1019, 990, 991, 12, 13, 268, 269, 524, 525, 780, 781, 46, 47, 28,
  94   29, 284, 285, 540, 541, 796, 797, 62, 63, 44, 45, 300, 301, 556, 557, 812,
  95   813, 302, 303, 60, 61, 316, 317, 572, 573, 828, 829, 318, 319, 76, 77,
  96   332, 333, 588, 589, 844, 845, 558, 559, 92, 93, 348, 349, 604, 605, 860,
  97   861, 574, 575, 108, 109, 364, 365, 620, 621, 876, 877, 814, 815, 124, 125,
  98   380, 381, 636, 637, 892, 893, 830, 831, 14, 15, 270, 271, 526, 527, 782,
  99   783, 110, 111, 30, 31, 286, 287, 542, 543, 798, 799, 126, 127, 140, 141,
 100   396, 397, 652, 653, 908, 909, 174, 175, 156, 157, 412, 413, 668, 669, 924,
 101   925, 190, 191, 172, 173, 428, 429, 684, 685, 940, 941, 430, 431, 188, 189,
 102   444, 445, 700, 701, 956, 957, 446, 447, 204, 205, 460, 461, 716, 717, 972,
 103   973, 686, 687, 220, 221, 476, 477, 732, 733, 988, 989, 702, 703, 236, 237,
 104   492, 493, 748, 749, 1004, 1005, 942, 943, 252, 253, 508, 509, 764, 765,
 105   1020, 1021, 958, 959, 142, 143, 398, 399, 654, 655, 910, 911, 238, 239, 158,
 106   159, 414, 415, 670, 671, 926, 927, 254, 255};
 107 #endif
 108
 109 /* construct a decimal64 NaN */
 110 static _Decimal64
 111 get_decimal64_nan (void)
 112 {
 113   union ieee_double_extract x;
 114   union ieee_double_decimal64 y;
 115
 116   x.s.exp = 1984; /* G[0]..G[4] = 11111: quiet NaN */
 117   y.d = x.d;
 118   return y.d64;
 119 }
 120
 121 /* construct the decimal64 Inf with given sign */
 122 static _Decimal64
 123 get_decimal64_inf (int negative)
 124 {
 125   union ieee_double_extract x;
 126   union ieee_double_decimal64 y;
 127
 128   x.s.sig = (negative) ? 1 : 0;
 129   x.s.exp = 1920; /* G[0]..G[4] = 11110: Inf */
 130   y.d = x.d;
 131   return y.d64;
 132 }
 133
 134 /* construct the decimal64 zero with given sign */
 135 static _Decimal64
 136 get_decimal64_zero (int negative)
 137 {
 138   union ieee_double_decimal64 y;
 139
 140   /* zero has the same representation in binary64 and decimal64 */
 141   y.d = negative ? DBL_NEG_ZERO : 0.0;
 142   return y.d64;
 143 }
 144
 145 /* construct the decimal64 smallest non-zero with given sign */
 146 static _Decimal64
 147 get_decimal64_min (int negative)
 148 {
 149   return negative ? - 1E-398dd : 1E-398dd;
 150 }
 151
 152 /* construct the decimal64 largest finite number with given sign */
 153 static _Decimal64
 154 get_decimal64_max (int negative)
 155 {
 156   return negative ? - DEC64_MAX : DEC64_MAX;
 157 }
 158
 159 /* one-to-one conversion:
 160    s is a decimal string representing a number x = m * 10^e which must be
 161    exactly representable in the decimal64 format, i.e.
 162    (a) the mantissa m has at most 16 decimal digits
 163    (b1) -383 <= e <= 384 with m integer multiple of 10^(-15), |m| < 10
 164    (b2) or -398 <= e <= 369 with m integer, |m| < 10^16.
 165    Assumes s is neither NaN nor +Inf nor -Inf.
 166 */
 167 static _Decimal64
 168 string_to_Decimal64 (char *s)
 169 {
 170   long int exp = 0;
 171   char m[17];
 172   long n = 0; /* mantissa length */
 173   char *endptr[1];
 174   union ieee_double_extract x;
 175   union ieee_double_decimal64 y;
 176 #ifdef DPD_FORMAT
 177   unsigned int G, d1, d2, d3, d4, d5;
 178 #endif
 179
 180   /* read sign */
 181   if (*s == '-')
 182     {
 183       x.s.sig = 1;
 184       s ++;
 185     }
 186   else
 187     x.s.sig = 0;
 188   /* read mantissa */
 189   while (ISDIGIT (*s))
 190     m[n++] = *s++;
 191   exp = n;
 192   if (*s == '.')
 193     {
 194       s ++;
 195       while (ISDIGIT (*s))
 196         m[n++] = *s++;
 197     }
 198   /* we have exp digits before decimal point, and a total of n digits */
 199   exp -= n; /* we will consider an integer mantissa */
 200   MPFR_ASSERTN(n <= 16);
 201   if (*s == 'E' || *s == 'e')
 202     exp += strtol (s + 1, endptr, 10);
 203   else
 204     *endptr = s;
 205   MPFR_ASSERTN(**endptr == '\0');
 206   MPFR_ASSERTN(-398 <= exp && exp <= (long) (385 - n));
 207   while (n < 16)
 208     {
 209       m[n++] = '0';
 210       exp --;
 211     }
 212   /* now n=16 and -398 <= exp <= 369 */
 213   m[n] = '\0';
 214
 215   /* compute biased exponent */
 216   exp += 398;
 217
 218   MPFR_ASSERTN(exp >= -15);
 219   if (exp < 0)
 220     {
 221       int i;
 222       n = -exp;
 223       /* check the last n digits of the mantissa are zero */
 224       for (i = 1; i <= n; i++)
 225         MPFR_ASSERTN(m[16 - n] == '0');
 226       /* shift the first (16-n) digits to the right */
 227       for (i = 16 - n - 1; i >= 0; i--)
 228         m[i + n] = m[i];
 229       /* zero the first n digits */
 230       for (i = 0; i < n; i ++)
 231         m[i] = '0';
 232       exp = 0;
 233     }
 234
 235   /* now convert to DPD or BID */
 236 #ifdef DPD_FORMAT
 237 #define CH(d) (d - '0')
 238   if (m[0] >= '8')
 239     G = (3 << 11) | ((exp & 768) << 1) | ((CH(m[0]) & 1) << 8);
 240   else
 241     G = ((exp & 768) << 3) | (CH(m[0]) << 8);
 242   /* now the most 5 significant bits of G are filled */
 243   G |= exp & 255;
 244   d1 = T[100 * CH(m[1]) + 10 * CH(m[2]) + CH(m[3])]; /* 10-bit encoding */
 245   d2 = T[100 * CH(m[4]) + 10 * CH(m[5]) + CH(m[6])]; /* 10-bit encoding */
 246   d3 = T[100 * CH(m[7]) + 10 * CH(m[8]) + CH(m[9])]; /* 10-bit encoding */
 247   d4 = T[100 * CH(m[10]) + 10 * CH(m[11]) + CH(m[12])]; /* 10-bit encoding */
 248   d5 = T[100 * CH(m[13]) + 10 * CH(m[14]) + CH(m[15])]; /* 10-bit encoding */
 249   x.s.exp = G >> 2;
 250   x.s.manh = ((G & 3) << 18) | (d1 << 8) | (d2 >> 2);
 251   x.s.manl = (d2 & 3) << 30;
 252   x.s.manl |= (d3 << 20) | (d4 << 10) | d5;
 253 #else /* BID format */
 254   {
 255     mp_size_t rn;
 256     mp_limb_t rp[2];
 257     int case_i = strcmp (m, "9007199254740992") < 0;
 258
 259     for (n = 0; n < 16; n++)
 260       m[n] -= '0';
 261     rn = mpn_set_str (rp, (unsigned char *) m, 16, 10);
 262     if (rn == 1)
 263       rp[1] = 0;
 264 #if GMP_NUMB_BITS > 32
 265     rp[1] = rp[1] << (GMP_NUMB_BITS - 32);
 266     rp[1] |= rp[0] >> 32;
 267     rp[0] &= 4294967295UL;
 268 #endif
 269     if (case_i)
 270       {  /* s < 2^53: case i) */
 271         x.s.exp = exp << 1;
 272         x.s.manl = rp[0];           /* 32 bits */
 273         x.s.manh = rp[1] & 1048575; /* 20 low bits */
 274         x.s.exp |= rp[1] >> 20;     /* 1 bit */
 275       }
 276     else /* s >= 2^53: case ii) */
 277       {
 278         x.s.exp = 1536 | (exp >> 1);
 279         x.s.manl = rp[0];
 280         x.s.manh = (rp[1] ^ 2097152) | ((exp & 1) << 19);
 281       }
 282   }
 283 #endif /* DPD_FORMAT */
 284   y.d = x.d;
 285   return y.d64;
 286 }
 287
 288 _Decimal64
 289 mpfr_get_decimal64 (mpfr_srcptr src, mpfr_rnd_t rnd_mode)
 290 {
 291   int negative;
 292   mpfr_exp_t e;
 293
 294   /* the encoding of NaN, Inf, zero is the same under DPD or BID */
 295   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
 296     {
 297       if (MPFR_IS_NAN (src))
 298         return get_decimal64_nan ();
 299
 300       negative = MPFR_IS_NEG (src);
 301
 302       if (MPFR_IS_INF (src))
 303         return get_decimal64_inf (negative);
 304
 305       MPFR_ASSERTD (MPFR_IS_ZERO(src));
 306       return get_decimal64_zero (negative);
 307     }
 308
 309   e = MPFR_GET_EXP (src);
 310   negative = MPFR_IS_NEG (src);
 311
 312   if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDA))
 313     rnd_mode = negative ? MPFR_RNDD : MPFR_RNDU;
 314
 315   /* the smallest decimal64 number is 10^(-398),
 316      with 2^(-1323) < 10^(-398) < 2^(-1322) */
 317   if (MPFR_UNLIKELY (e < -1323)) /* src <= 2^(-1324) < 1/2*10^(-398) */
 318     {
 319       if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDN
 320           || (rnd_mode == MPFR_RNDD && negative == 0)
 321           || (rnd_mode == MPFR_RNDU && negative != 0))
 322         return get_decimal64_zero (negative);
 323       else /* return the smallest non-zero number */
 324         return get_decimal64_min (negative);
 325     }
 326   /* the largest decimal64 number is just below 10^(385) < 2^1279 */
 327   else if (MPFR_UNLIKELY (e > 1279)) /* then src >= 2^1279 */
 328     {
 329       if (rnd_mode == MPFR_RNDZ
 330           || (rnd_mode == MPFR_RNDU && negative != 0)
 331           || (rnd_mode == MPFR_RNDD && negative == 0))
 332         return get_decimal64_max (negative);
 333       else
 334         return get_decimal64_inf (negative);
 335     }
 336   else
 337     {
 338       /* we need to store the sign (1), the mantissa (16), and the terminating
 339          character, thus we need at least 18 characters in s */
 340       char s[23];
 341       mpfr_get_str (s, &e, 10, 16, src, rnd_mode);
 342       /* the smallest normal number is 1.000...000E-383,
 343          which corresponds to s=[0.]1000...000 and e=-382 */
 344       if (e < -382)
 345         {
 346           /* the smallest subnormal number is 0.000...001E-383 = 1E-398,
 347              which corresponds to s=[0.]1000...000 and e=-397 */
 348           if (e < -397)
 349             {
 350               if (rnd_mode == MPFR_RNDN && e == -398)
 351                 {
 352                   /* If 0.5E-398 < |src| < 1E-398 (smallest subnormal),
 353                      src should round to +/- 1E-398 in MPFR_RNDN. */
 354                   mpfr_get_str (s, &e, 10, 1, src, MPFR_RNDA);
 355                   return e == -398 && s[negative] <= '5' ?
 356                     get_decimal64_zero (negative) :
 357                     get_decimal64_min (negative);
 358                 }
 359               if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDN
 360                   || (rnd_mode == MPFR_RNDD && negative == 0)
 361                   || (rnd_mode == MPFR_RNDU && negative != 0))
 362                 return get_decimal64_zero (negative);
 363               else /* return the smallest non-zero number */
 364                 return get_decimal64_min (negative);
 365             }
 366           else
 367             {
 368               mpfr_exp_t e2;
 369               long digits = 16 - (-382 - e);
 370               /* if e = -397 then 16 - (-382 - e) = 1 */
 371               mpfr_get_str (s, &e2, 10, digits, src, rnd_mode);
 372               /* Warning: we can have e2 = e + 1 here, when rounding to
 373                  nearest or away from zero. */
 374               s[negative + digits] = 'E';
 375               sprintf (s + negative + digits + 1, "%ld",
 376                        (long int)e2 - digits);
 377               return string_to_Decimal64 (s);
 378             }
 379         }
 380       /* the largest number is 9.999...999E+384,
 381          which corresponds to s=[0.]9999...999 and e=385 */
 382       else if (e > 385)
 383         {
 384           if (rnd_mode == MPFR_RNDZ
 385               || (rnd_mode == MPFR_RNDU && negative != 0)
 386               || (rnd_mode == MPFR_RNDD && negative == 0))
 387             return get_decimal64_max (negative);
 388           else
 389             return get_decimal64_inf (negative);
 390         }
 391       else /* -382 <= e <= 385 */
 392         {
 393           s[16 + negative] = 'E';
 394           sprintf (s + 17 + negative, "%ld", (long int)e - 16);
 395           return string_to_Decimal64 (s);
 396         }
 397     }
 398 }
 399
 400 #endif /* MPFR_WANT_DECIMAL_FLOATS */