contrib/mpfr/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 Free Software Foundation, Inc.
   9 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  25 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  26 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 BITS_PER_MP_LIMB > 32
 273     rp[1] = rp[1] << (BITS_PER_MP_LIMB - 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, mp_rnd_t rnd_mode)
 298 {
 299   int negative;
 300   mp_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   /* the smallest decimal64 number is 10^(-398),
 321      with 2^(-1323) < 10^(-398) < 2^(-1322) */
 322   if (MPFR_UNLIKELY (e < -1323)) /* src <= 2^(-1324) < 1/2*10^(-398) */
 323     {
 324       if (rnd_mode == GMP_RNDZ || rnd_mode == GMP_RNDN
 325           || (rnd_mode == GMP_RNDD && negative == 0)
 326           || (rnd_mode == GMP_RNDU && negative != 0))
 327         return get_decimal64_zero (negative);
 328       else /* return the smallest non-zero number */
 329         return get_decimal64_min (negative);
 330     }
 331   /* the largest decimal64 number is just below 10^(385) < 2^1279 */
 332   else if (MPFR_UNLIKELY (e > 1279)) /* then src >= 2^1279 */
 333     {
 334       if (GMP_RNDZ || (rnd_mode == GMP_RNDU && negative != 0)
 335           || (rnd_mode == GMP_RNDD && negative == 0))
 336         return get_decimal64_max (negative);
 337       else
 338         return get_decimal64_inf (negative);
 339     }
 340   else
 341     {
 342       /* we need to store the sign (1), the mantissa (16), and the terminating
 343          character, thus we need at least 18 characters in s */
 344       char s[23];
 345       mpfr_get_str (s, &e, 10, 16, src, rnd_mode);
 346       /* the smallest normal number is 1.000...000E-383,
 347          which corresponds to s=[0.]1000...000 and e=-382 */
 348       if (e < -382)
 349         {
 350           /* the smallest subnormal number is 0.000...001E-383 = 1E-398,
 351              which corresponds to s=[0.]1000...000 and e=-397 */
 352           if (e < -397)
 353             {
 354               if (rnd_mode == GMP_RNDZ || rnd_mode == GMP_RNDN
 355                   || (rnd_mode == GMP_RNDD && negative == 0)
 356                   || (rnd_mode == GMP_RNDU && negative != 0))
 357                 return get_decimal64_zero (negative);
 358               else /* return the smallest non-zero number */
 359                 return get_decimal64_min (negative);
 360             }
 361           else
 362             {
 363               mp_exp_t e2;
 364               long digits = 16 - (-382 - e);
 365               /* if e = -397 then 16 - (-382 - e) = 1 */
 366               mpfr_get_str (s, &e2, 10, digits, src, rnd_mode);
 367               /* Warning: we can have e2 = e + 1 here, when rounding to
 368                  nearest or away from zero. */
 369               s[negative + digits] = 'E';
 370               sprintf (s + negative + digits + 1, "%ld",
 371                        (long int)e2 - digits);
 372               return string_to_Decimal64 (s);
 373             }
 374         }
 375       /* the largest number is 9.999...999E+384,
 376          which corresponds to s=[0.]9999...999 and e=385 */
 377       else if (e > 385)
 378         {
 379           if (GMP_RNDZ || (rnd_mode == GMP_RNDU && negative != 0)
 380               || (rnd_mode == GMP_RNDD && negative == 0))
 381             return get_decimal64_max (negative);
 382           else
 383             return get_decimal64_inf (negative);
 384         }
 385       else /* -382 <= e <= 385 */
 386         {
 387           s[16 + negative] = 'E';
 388           sprintf (s + 17 + negative, "%ld", (long int)e - 16);
 389           return string_to_Decimal64 (s);
 390         }
 391     }
 392 }
 393
 394 #endif /* MPFR_WANT_DECIMAL_FLOATS */