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