Home | History | Annotate | Line # | Download | only in src 
      1 /* mpfr_get_decimal128 -- convert a multiple precision floating-point number
      2                           to an IEEE 754-2008 decimal128 float
      3 
      4 See https://gcc.gnu.org/legacy-ml/gcc/2006-06/msg00691.html,
      5 https://gcc.gnu.org/onlinedocs/gcc/Decimal-Float.html,
      6 and TR 24732 <https://www.open-std.org/jtc1/sc22/wg14/www/projects#24732>.
      7 
      8 Copyright 2006-2023 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 https://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 /* Warning! Do not use any conversion between binary and decimal types,
     29  * otherwise GCC will generate from 2 to 3 MB of code (depending on the
     30  * GCC version) in the MPFR shared library when the _Decimal128 format
     31  * is BID (e.g. on x86).
     32  *   https://gcc.gnu.org/bugzilla/show_bug.cgi?id=96173
     33  *   https://gforge.inria.fr/tracker/index.php?func=detail&aid=21849&group_id=136&atid=619
     34  *
     35  * FIXME: Try to save even more space in the MPFR library by avoiding
     36  * _Decimal128 operations entirely. These operations now appear only in
     37  * string_to_Decimal128(). In the case where the _Decimal128 format is
     38  * recognized as BID, this function should be reimplemented directly by
     39  * using the specification of the encoding of this format, as already
     40  * done for _Decimal64 (see string_to_Decimal64 in get_d64.c).
     41  * Or use strtod128 when available, making sure that the string is
     42  * locale-independent? (Should one optionally use libdfp for that?)
     43  */
     44 
     45 #include "mpfr-impl.h"
     46 #include "ieee_floats.h"
     47 
     48 #define ISDIGIT(c) ('0' <= c && c <= '9')
     49 
     50 #ifdef MPFR_WANT_DECIMAL_FLOATS
     51 
     52 #ifndef DEC128_MAX
     53 # define DEC128_MAX 9.999999999999999999999999999999999E6144dl
     54 #endif
     55 
     56 /* construct a decimal128 NaN */
     57 static _Decimal128
     58 get_decimal128_nan (void)
     59 {
     60   return 0.dl / 0.dl;
     61 }
     62 
     63 /* construct the decimal128 Inf with given sign */
     64 static _Decimal128
     65 get_decimal128_inf (int negative)
     66 {
     67   return negative ? - 1.dl / 0.dl : 1.dl / 0.dl;
     68 }
     69 
     70 /* construct the decimal128 zero with given sign */
     71 static _Decimal128
     72 get_decimal128_zero (int negative)
     73 {
     74   return negative ? - 0.dl : 0.dl;
     75 }
     76 
     77 /* construct the decimal128 smallest non-zero with given sign:
     78    it is 10^emin * 10^(1-p). Since emax = 6144, emin = 1-emax = -6143,
     79    and p = 34, we get 10^(-6176) */
     80 static _Decimal128
     81 get_decimal128_min (int negative)
     82 {
     83   return negative ? - 1E-6176dl : 1E-6176dl;
     84 }
     85 
     86 /* construct the decimal128 largest finite number with given sign */
     87 static _Decimal128
     88 get_decimal128_max (int negative)
     89 {
     90   return negative ? - DEC128_MAX : DEC128_MAX;
     91 }
     92 
     93 /* one-to-one conversion:
     94    s is a decimal string representing a number x = m * 10^e which must be
     95    exactly representable in the decimal128 format, i.e.
     96    (a) the mantissa m has at most 34 decimal digits
     97    (b1) -6143 <= e <= 6144 with m integer multiple of 10^(-33), |m| < 10
     98    (b2) or -6176 <= e <= 6111 with m integer, |m| < 10^34.
     99    Assumes s is neither NaN nor +Inf nor -Inf.
    100    s = [-][0-9]+E[-][0-9]+
    101 
    102    The decimal128 format (cf table 3.6 of IEEE 754-2008) has the following
    103    parameters:
    104    * k = 128 (number of bits of storage)
    105    * p = 34 (precision in digits)
    106    * emax = 6144
    107    * bias = E-q = 6176
    108    * sign bit has 1 bit
    109    * w+5 = 17 bits (combination field width)
    110    * t = 110 bits (trailing significand width)
    111    We have k = 1 + 5 + w + t = 128.
    112 */
    113 static _Decimal128
    114 string_to_Decimal128 (char *s) /* portable version */
    115 {
    116   long int exp = 0;
    117   char m[35];
    118   long n = 0; /* mantissa length */
    119   char *endptr[1];
    120   _Decimal128 x = 0;
    121   int sign = 0;
    122 
    123   /* read sign */
    124   if (*s == '-')
    125     {
    126       sign = 1;
    127       s ++;
    128     }
    129   /* read mantissa */
    130   while (ISDIGIT (*s))
    131     m[n++] = *s++;
    132 
    133   /* as constructed in mpfr_get_decimal128, s cannot have any '.' separator */
    134 
    135   /* we will consider an integer mantissa m*10^exp */
    136   MPFR_ASSERTN(n <= 34);
    137   /* s always has an exponent separator 'E' */
    138   MPFR_ASSERTN(*s == 'E');
    139   exp = strtol (s + 1, endptr, 10);
    140   MPFR_ASSERTN(**endptr == '\0');
    141   MPFR_ASSERTN(-6176 <= exp && exp <= (long) (6145 - n));
    142   while (n < 34)
    143     {
    144       m[n++] = '0';
    145       exp --;
    146     }
    147   /* now n=34 and -6176 <= exp <= 6111, cf (b2) */
    148   m[n] = '\0';
    149 
    150   /* the number to convert is m[] * 10^exp where the mantissa is a 34-digit
    151      integer */
    152 
    153   /* compute biased exponent */
    154   exp += 6176;
    155 
    156   MPFR_ASSERTN(exp >= -33);
    157   if (exp < 0)
    158     {
    159       int i;
    160       n = -exp;
    161       /* check the last n digits of the mantissa are zero */
    162       for (i = 1; i <= n; i++)
    163         MPFR_ASSERTN(m[34 - n] == '0');
    164       /* shift the first (34-n) digits to the right */
    165       for (i = 34 - n - 1; i >= 0; i--)
    166         m[i + n] = m[i];
    167       /* zero the first n digits */
    168       for (i = 0; i < n; i ++)
    169         m[i] = '0';
    170       exp = 0;
    171     }
    172 
    173   /* the number to convert is m[] * 10^(exp-6176) */
    174   exp -= 6176;
    175 
    176   for (n = 0; n < 34; n++)
    177     x = (_Decimal128) 10 * x + (_Decimal128) (m[n] - '0');
    178 
    179   /* multiply by 10^exp */
    180   if (exp > 0)
    181     {
    182       _Decimal128 ten = 10;
    183       _Decimal128 ten2 = ten * ten;
    184       _Decimal128 ten4 = ten2 * ten2;
    185       _Decimal128 ten8 = ten4 * ten4;
    186       _Decimal128 ten16 = ten8 * ten8;
    187       _Decimal128 ten32 = ten16 * ten16;
    188       _Decimal128 ten64 = ten32 * ten32;
    189       _Decimal128 ten128 = ten64 * ten64;
    190       _Decimal128 ten256 = ten128 * ten128;
    191       _Decimal128 ten512 = ten256 * ten256;
    192       _Decimal128 ten1024 = ten512 * ten512;
    193       _Decimal128 ten2048 = ten1024 * ten1024;
    194       _Decimal128 ten4096 = ten2048 * ten2048;
    195 
    196       if (exp >= 4096)
    197         {
    198           x *= ten4096;
    199           exp -= 4096;
    200         }
    201       if (exp >= 2048)
    202         {
    203           x *= ten2048;
    204           exp -= 2048;
    205         }
    206       if (exp >= 1024)
    207         {
    208           x *= ten1024;
    209           exp -= 1024;
    210         }
    211       if (exp >= 512)
    212         {
    213           x *= ten512;
    214           exp -= 512;
    215         }
    216       if (exp >= 256)
    217         {
    218           x *= ten256;
    219           exp -= 256;
    220         }
    221       if (exp >= 128)
    222         {
    223           x *= ten128;
    224           exp -= 128;
    225         }
    226       if (exp >= 64)
    227         {
    228           x *= ten64;
    229           exp -= 64;
    230         }
    231       if (exp >= 32)
    232         {
    233           x *= ten32;
    234           exp -= 32;
    235         }
    236       if (exp >= 16)
    237         {
    238           x *= ten16;
    239           exp -= 16;
    240         }
    241       if (exp >= 8)
    242         {
    243           x *= ten8;
    244           exp -= 8;
    245         }
    246       if (exp >= 4)
    247         {
    248           x *= ten4;
    249           exp -= 4;
    250         }
    251       if (exp >= 2)
    252         {
    253           x *= ten2;
    254           exp -= 2;
    255         }
    256       if (exp >= 1)
    257         {
    258           x *= ten;
    259           exp -= 1;
    260         }
    261     }
    262   else if (exp < 0)
    263     {
    264       _Decimal128 ten = 10;
    265       _Decimal128 ten2 = ten * ten;
    266       _Decimal128 ten4 = ten2 * ten2;
    267       _Decimal128 ten8 = ten4 * ten4;
    268       _Decimal128 ten16 = ten8 * ten8;
    269       _Decimal128 ten32 = ten16 * ten16;
    270       _Decimal128 ten64 = ten32 * ten32;
    271       _Decimal128 ten128 = ten64 * ten64;
    272       _Decimal128 ten256 = ten128 * ten128;
    273       _Decimal128 ten512 = ten256 * ten256;
    274       _Decimal128 ten1024 = ten512 * ten512;
    275       _Decimal128 ten2048 = ten1024 * ten1024;
    276       _Decimal128 ten4096 = ten2048 * ten2048;
    277 
    278       if (exp <= -4096)
    279         {
    280           x /= ten4096;
    281           exp += 4096;
    282         }
    283       if (exp <= -2048)
    284         {
    285           x /= ten2048;
    286           exp += 2048;
    287         }
    288       if (exp <= -1024)
    289         {
    290           x /= ten1024;
    291           exp += 1024;
    292         }
    293       if (exp <= -512)
    294         {
    295           x /= ten512;
    296           exp += 512;
    297         }
    298       if (exp <= -256)
    299         {
    300           x /= ten256;
    301           exp += 256;
    302         }
    303       if (exp <= -128)
    304         {
    305           x /= ten128;
    306           exp += 128;
    307         }
    308       if (exp <= -64)
    309         {
    310           x /= ten64;
    311           exp += 64;
    312         }
    313       if (exp <= -32)
    314         {
    315           x /= ten32;
    316           exp += 32;
    317         }
    318       if (exp <= -16)
    319         {
    320           x /= ten16;
    321           exp += 16;
    322         }
    323       if (exp <= -8)
    324         {
    325           x /= ten8;
    326           exp += 8;
    327         }
    328       if (exp <= -4)
    329         {
    330           x /= ten4;
    331           exp += 4;
    332         }
    333       if (exp <= -2)
    334         {
    335           x /= ten2;
    336           exp += 2;
    337         }
    338       if (exp <= -1)
    339         {
    340           x /= ten;
    341           exp += 1;
    342         }
    343     }
    344 
    345   if (sign)
    346     x = -x;
    347 
    348   return x;
    349 }
    350 
    351 _Decimal128
    352 mpfr_get_decimal128 (mpfr_srcptr src, mpfr_rnd_t rnd_mode)
    353 {
    354   int negative;
    355   mpfr_exp_t e;
    356 
    357   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
    358     {
    359       if (MPFR_IS_NAN (src))
    360         {
    361           /* we don't propagate the sign bit */
    362           return get_decimal128_nan ();
    363         }
    364 
    365       negative = MPFR_IS_NEG (src);
    366 
    367       if (MPFR_IS_INF (src))
    368         return get_decimal128_inf (negative);
    369 
    370       MPFR_ASSERTD (MPFR_IS_ZERO(src));
    371       return get_decimal128_zero (negative);
    372     }
    373 
    374   e = MPFR_GET_EXP (src);
    375   negative = MPFR_IS_NEG (src);
    376 
    377   MPFR_UPDATE2_RND_MODE (rnd_mode, MPFR_SIGN (src));
    378 
    379   /* now rnd_mode is RNDN, RNDF, RNDA or RNDZ */
    380 
    381   /* the smallest decimal128 number is 10^(-6176),
    382      with 2^(-20517) < 10^(-6176) < 2^(-20516) */
    383   if (MPFR_UNLIKELY (e < -20517)) /* src <= 2^(-20518) < 1/2*10^(-6176) */
    384     {
    385       if (rnd_mode != MPFR_RNDA)
    386         return get_decimal128_zero (negative);
    387       else /* RNDA: return the smallest non-zero number */
    388         return get_decimal128_min (negative);
    389     }
    390   /* the largest decimal128 number is just below 10^6145 < 2^20414 */
    391   else if (MPFR_UNLIKELY (e > 20414)) /* then src >= 2^20414 */
    392     {
    393       if (rnd_mode == MPFR_RNDZ)
    394         return get_decimal128_max (negative);
    395       else /* RNDN, RNDA, RNDF: round away */
    396         return get_decimal128_inf (negative);
    397     }
    398   else
    399     {
    400       /* we need to store the sign (1 character), the significand (at most 34
    401          characters), the exponent part (at most 6 characters for "E-6176"),
    402          and the terminating null character, thus we need at least 42
    403          characters in s. */
    404       char s[42];
    405       mpfr_get_str (s, &e, 10, 34, src, rnd_mode);
    406       /* the smallest normal number is 1.000...000E-6143,
    407          which corresponds to s=[0.]1000...000 and e=-6142 */
    408       if (e < -6142)
    409         {
    410           /* the smallest subnormal number is 0.000...001E-6143 = 1E-6176,
    411              which corresponds to s=[0.]1000...000 and e=-6175 */
    412           if (e < -6175)
    413             {
    414               if (rnd_mode == MPFR_RNDN && e == -6176)
    415                 {
    416                   /* If 0.5E-6176 < |src| < 1E-6176 (smallest subnormal),
    417                      src should round to +/- 1E-6176 in MPFR_RNDN. */
    418                   mpfr_get_str (s, &e, 10, 1, src, MPFR_RNDA);
    419                   return e == -6176 && s[negative] <= '5' ?
    420                     get_decimal128_zero (negative) :
    421                     get_decimal128_min (negative);
    422                 }
    423               if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDN)
    424                 return get_decimal128_zero (negative);
    425               else /* RNDA or RNDF: return the smallest non-zero number */
    426                 return get_decimal128_min (negative);
    427             }
    428           else
    429             {
    430               mpfr_exp_t e2;
    431               long digits = 34 - (-6142 - e);
    432               /* if e = -6175 then 34 - (-6142 - e) = 1 */
    433               mpfr_get_str (s, &e2, 10, digits, src, rnd_mode);
    434               /* Warning: we can have e2 = e + 1 here, when rounding to
    435                  nearest or away from zero. */
    436               s[negative + digits] = 'E';
    437               sprintf (s + negative + digits + 1, "%ld",
    438                        (long int)e2 - digits);
    439               return string_to_Decimal128 (s);
    440             }
    441         }
    442       /* the largest number is 9.999...999E+6144,
    443          which corresponds to s=[0.]9999...999 and e=6145 */
    444       else if (e > 6145)
    445         {
    446           if (rnd_mode == MPFR_RNDZ)
    447             return get_decimal128_max (negative);
    448           else /* RNDN, RNDA, RNDF: round away */
    449             return get_decimal128_inf (negative);
    450         }
    451       else /* -6142 <= e <= 6145 */
    452         {
    453           s[34 + negative] = 'E';
    454           sprintf (s + 35 + negative, "%ld", (long int) e - 34);
    455           return string_to_Decimal128 (s);
    456         }
    457     }
    458 }
    459 
    460 #endif /* MPFR_WANT_DECIMAL_FLOATS */
    461