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