lib/strtod.c

   1 /* Copyright (C) 1991-1992, 1997, 1999, 2003, 2006, 2008-2024 Free Software
   2    Foundation, Inc.
   3
   4    This file is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU Lesser General Public License as
   6    published by the Free Software Foundation, either version 3 of the
   7    License, or (at your option) any later version.
   8
   9    This file is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU Lesser General Public License for more details.
  13
  14    You should have received a copy of the GNU Lesser General Public License
  15    along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
  16
  17 #if ! (defined USE_FLOAT || defined USE_LONG_DOUBLE)
  18 # include <config.h>
  19 #endif
  20
  21 /* Specification.  */
  22 #include <stdlib.h>
  23
  24 #include <ctype.h>      /* isspace() */
  25 #include <errno.h>
  26 #include <float.h>      /* {FLT,DBL,LDBL}_{MIN,MAX} */
  27 #include <limits.h>     /* LONG_{MIN,MAX} */
  28 #include <locale.h>     /* localeconv() */
  29 #include <math.h>       /* NAN */
  30 #include <stdio.h>      /* sprintf() */
  31 #include <string.h>     /* strdup() */
  32 #if HAVE_NL_LANGINFO
  33 # include <langinfo.h>
  34 #endif
  35
  36 #include "c-ctype.h"
  37
  38 #undef MIN
  39 #undef MAX
  40 #if defined USE_FLOAT
  41 # define STRTOD strtof
  42 # define LDEXP ldexpf
  43 # if STRTOF_HAS_UNDERFLOW_BUG
  44    /* strtof would not set errno=ERANGE upon flush-to-zero underflow.  */
  45 #  define HAVE_UNDERLYING_STRTOD 0
  46 # else
  47 #  define HAVE_UNDERLYING_STRTOD HAVE_STRTOF
  48 # endif
  49 # define HAS_GRADUAL_UNDERFLOW_PROBLEM STRTOF_HAS_GRADUAL_UNDERFLOW_PROBLEM
  50 # define DOUBLE float
  51 # define MIN FLT_MIN
  52 # define MAX FLT_MAX
  53 # define L_(literal) literal##f
  54 # if HAVE_LDEXPF_IN_LIBC
  55 #  define USE_LDEXP 1
  56 # else
  57 #  define USE_LDEXP 0
  58 # endif
  59 #elif defined USE_LONG_DOUBLE
  60 # define STRTOD strtold
  61 # define LDEXP ldexpl
  62 # if defined __hpux && defined __hppa
  63    /* We cannot call strtold on HP-UX/hppa, because its return type is a struct,
  64       not a 'long double'.  */
  65 #  define HAVE_UNDERLYING_STRTOD 0
  66 # elif STRTOLD_HAS_UNDERFLOW_BUG
  67    /* strtold would not set errno=ERANGE upon flush-to-zero underflow.  */
  68 #  define HAVE_UNDERLYING_STRTOD 0
  69 # elif defined __MINGW32__ && __MINGW64_VERSION_MAJOR < 10
  70    /* strtold is broken in mingw versions before 10.0:
  71       - Up to mingw 5.0.x, it leaks memory at every invocation.
  72       - Up to mingw 9.0.x, it allocates an unbounded amount of stack.
  73       See <https://github.com/mingw-w64/mingw-w64/commit/450309b97b2e839ea02887dfaf0f1d10fb5d40cc>
  74       and <https://github.com/mingw-w64/mingw-w64/commit/73806c0709b7e6c0f6587f11a955743670e85470>.  */
  75 #  define HAVE_UNDERLYING_STRTOD 0
  76 # elif defined __HAIKU__
  77    /* Haiku's strtold maps denormalized numbers to zero.
  78       <https://dev.haiku-os.org/ticket/19040>  */
  79 #  define HAVE_UNDERLYING_STRTOD 0
  80 # else
  81 #  define HAVE_UNDERLYING_STRTOD HAVE_STRTOLD
  82 # endif
  83 # define HAS_GRADUAL_UNDERFLOW_PROBLEM STRTOLD_HAS_GRADUAL_UNDERFLOW_PROBLEM
  84 # define DOUBLE long double
  85 # define MIN LDBL_MIN
  86 # define MAX LDBL_MAX
  87 # define L_(literal) literal##L
  88 # if HAVE_LDEXPL_IN_LIBC
  89 #  define USE_LDEXP 1
  90 # else
  91 #  define USE_LDEXP 0
  92 # endif
  93 #else
  94 # define STRTOD strtod
  95 # define LDEXP ldexp
  96 # if STRTOD_HAS_UNDERFLOW_BUG
  97    /* strtod would not set errno=ERANGE upon flush-to-zero underflow.  */
  98 #  define HAVE_UNDERLYING_STRTOD 0
  99 # else
 100 #  define HAVE_UNDERLYING_STRTOD 1
 101 # endif
 102 # define HAS_GRADUAL_UNDERFLOW_PROBLEM STRTOD_HAS_GRADUAL_UNDERFLOW_PROBLEM
 103 # define DOUBLE double
 104 # define MIN DBL_MIN
 105 # define MAX DBL_MAX
 106 # define L_(literal) literal
 107 # if HAVE_LDEXP_IN_LIBC
 108 #  define USE_LDEXP 1
 109 # else
 110 #  define USE_LDEXP 0
 111 # endif
 112 #endif
 113
 114 /* Return true if C is a space in the current locale, avoiding
 115    problems with signed char and isspace.  */
 116 static bool
 117 locale_isspace (char c)
 118 {
 119   unsigned char uc = c;
 120   return isspace (uc) != 0;
 121 }
 122
 123 /* Determine the decimal-point character according to the current locale.  */
 124 static char
 125 decimal_point_char (void)
 126 {
 127   const char *point;
 128   /* Determine it in a multithread-safe way.  We know nl_langinfo is
 129      multithread-safe on glibc systems and Mac OS X systems, but is not required
 130      to be multithread-safe by POSIX.  sprintf(), however, is multithread-safe.
 131      localeconv() is rarely multithread-safe.  */
 132 #if HAVE_NL_LANGINFO && (__GLIBC__ || defined __UCLIBC__ || (defined __APPLE__ && defined __MACH__))
 133   point = nl_langinfo (RADIXCHAR);
 134 #elif 1
 135   char pointbuf[5];
 136   sprintf (pointbuf, "%#.0f", 1.0);
 137   point = &pointbuf[1];
 138 #else
 139   point = localeconv () -> decimal_point;
 140 #endif
 141   /* The decimal point is always a single byte: either '.' or ','.  */
 142   return (point[0] != '\0' ? point[0] : '.');
 143 }
 144
 145 #if !USE_LDEXP
 146  #undef LDEXP
 147  #define LDEXP dummy_ldexp
 148  /* A dummy definition that will never be invoked.  */
 149  static DOUBLE LDEXP (_GL_UNUSED DOUBLE x, _GL_UNUSED int exponent)
 150  {
 151    abort ();
 152    return L_(0.0);
 153  }
 154 #endif
 155
 156 /* Return X * BASE**EXPONENT.  Return an extreme value and set errno
 157    to ERANGE if underflow or overflow occurs.  */
 158 static DOUBLE
 159 scale_radix_exp (DOUBLE x, int radix, long int exponent)
 160 {
 161   /* If RADIX == 10, this code is neither precise nor fast; it is
 162      merely a straightforward and relatively portable approximation.
 163      If N == 2, this code is precise on a radix-2 implementation,
 164      albeit perhaps not fast if ldexp is not in libc.  */
 165
 166   long int e = exponent;
 167
 168   if (USE_LDEXP && radix == 2)
 169     return LDEXP (x, e < INT_MIN ? INT_MIN : INT_MAX < e ? INT_MAX : e);
 170   else
 171     {
 172       DOUBLE r = x;
 173
 174       if (r != 0)
 175         {
 176           if (e < 0)
 177             {
 178               for (;;)
 179                 {
 180                   if (e++ == 0)
 181                     {
 182                       if (r < MIN && r > -MIN)
 183                         /* Gradual underflow, resulting in a denormalized
 184                            number.  */
 185                         errno = ERANGE;
 186                       break;
 187                     }
 188                   r /= radix;
 189                   if (r == 0)
 190                     {
 191                       /* Flush-to-zero underflow.  */
 192                       errno = ERANGE;
 193                       break;
 194                     }
 195                 }
 196             }
 197           else
 198             {
 199               while (e-- != 0)
 200                 {
 201                   if (r < -MAX / radix)
 202                     {
 203                       errno = ERANGE;
 204                       return -HUGE_VAL;
 205                     }
 206                   else if (MAX / radix < r)
 207                     {
 208                       errno = ERANGE;
 209                       return HUGE_VAL;
 210                     }
 211                   else
 212                     r *= radix;
 213                 }
 214             }
 215         }
 216
 217       return r;
 218     }
 219 }
 220
 221 /* Parse a number at NPTR; this is a bit like strtol (NPTR, ENDPTR)
 222    except there are no leading spaces or signs or "0x", and ENDPTR is
 223    nonnull.  The number uses a base BASE (either 10 or 16) fraction, a
 224    radix RADIX (either 10 or 2) exponent, and exponent character
 225    EXPCHAR.  BASE is RADIX**RADIX_MULTIPLIER.  */
 226 static DOUBLE
 227 parse_number (const char *nptr,
 228               int base, int radix, int radix_multiplier, char radixchar,
 229               char expchar,
 230               char **endptr)
 231 {
 232   const char *s = nptr;
 233   const char *digits_start;
 234   const char *digits_end;
 235   const char *radixchar_ptr;
 236   long int exponent;
 237   DOUBLE num;
 238
 239   /* First, determine the start and end of the digit sequence.  */
 240   digits_start = s;
 241   radixchar_ptr = NULL;
 242   for (;; ++s)
 243     {
 244       if (base == 16 ? c_isxdigit (*s) : c_isdigit (*s))
 245         ;
 246       else if (radixchar_ptr == NULL && *s == radixchar)
 247         {
 248           /* Record that we have found the decimal point.  */
 249           radixchar_ptr = s;
 250         }
 251       else
 252         /* Any other character terminates the digit sequence.  */
 253         break;
 254     }
 255   digits_end = s;
 256   /* Now radixchar_ptr == NULL or
 257      digits_start <= radixchar_ptr < digits_end.  */
 258
 259   if (false)
 260     { /* Unoptimized.  */
 261       exponent =
 262         (radixchar_ptr != NULL
 263          ? - (long int) (digits_end - radixchar_ptr - 1)
 264          : 0);
 265     }
 266   else
 267     { /* Remove trailing zero digits.  This reduces rounding errors for
 268          inputs such as 1.0000000000 or 10000000000e-10.  */
 269       while (digits_end > digits_start)
 270         {
 271           if (digits_end - 1 == radixchar_ptr || *(digits_end - 1) == '0')
 272             digits_end--;
 273           else
 274             break;
 275         }
 276       exponent =
 277         (radixchar_ptr != NULL
 278          ? (digits_end > radixchar_ptr
 279             ? - (long int) (digits_end - radixchar_ptr - 1)
 280             : (long int) (radixchar_ptr - digits_end))
 281          : (long int) (s - digits_end));
 282     }
 283
 284   /* Then, convert the digit sequence to a number.  */
 285   {
 286     const char *dp;
 287     num = 0;
 288     for (dp = digits_start; dp < digits_end; dp++)
 289       if (dp != radixchar_ptr)
 290         {
 291           int digit;
 292
 293           /* Make sure that multiplication by BASE will not overflow.  */
 294           if (!(num <= MAX / base))
 295             {
 296               /* The value of the digit and all subsequent digits don't matter,
 297                  since we have already gotten as many digits as can be
 298                  represented in a 'DOUBLE'.  This doesn't necessarily mean that
 299                  the result will overflow: The exponent may reduce it to within
 300                  range.  */
 301               exponent +=
 302                 (digits_end - dp)
 303                 - (radixchar_ptr >= dp && radixchar_ptr < digits_end ? 1 : 0);
 304               break;
 305             }
 306
 307           /* Eat the next digit.  */
 308           if (c_isdigit (*dp))
 309             digit = *dp - '0';
 310           else if (base == 16 && c_isxdigit (*dp))
 311             digit = c_tolower (*dp) - ('a' - 10);
 312           else
 313             abort ();
 314           num = num * base + digit;
 315         }
 316   }
 317
 318   exponent = exponent * radix_multiplier;
 319
 320   /* Finally, parse the exponent.  */
 321   if (c_tolower (*s) == expchar && ! locale_isspace (s[1]))
 322     {
 323       /* Add any given exponent to the implicit one.  */
 324       int saved_errno = errno;
 325       char *end;
 326       long int value = strtol (s + 1, &end, 10);
 327       errno = saved_errno;
 328
 329       if (s + 1 != end)
 330         {
 331           /* Skip past the exponent, and add in the implicit exponent,
 332              resulting in an extreme value on overflow.  */
 333           s = end;
 334           exponent =
 335             (exponent < 0
 336              ? (value < LONG_MIN - exponent ? LONG_MIN : exponent + value)
 337              : (LONG_MAX - exponent < value ? LONG_MAX : exponent + value));
 338         }
 339     }
 340
 341   *endptr = (char *) s;
 342   return scale_radix_exp (num, radix, exponent);
 343 }
 344
 345 /* HP cc on HP-UX 10.20 has a bug with the constant expression -0.0.
 346    ICC 10.0 has a bug when optimizing the expression -zero.
 347    The expression -MIN * MIN does not work when cross-compiling
 348    to PowerPC on Mac OS X 10.5.  */
 349 static DOUBLE
 350 minus_zero (void)
 351 {
 352 #if defined __hpux || defined __sgi || defined __ICC
 353   return -MIN * MIN;
 354 #else
 355   return -0.0;
 356 #endif
 357 }
 358
 359 /* Convert NPTR to a DOUBLE.  If ENDPTR is not NULL, a pointer to the
 360    character after the last one used in the number is put in *ENDPTR.  */
 361 DOUBLE
 362 STRTOD (const char *nptr, char **endptr)
 363 #if HAVE_UNDERLYING_STRTOD
 364 # if defined USE_FLOAT
 365 #  undef strtof
 366 # elif defined USE_LONG_DOUBLE
 367 #  undef strtold
 368 # else
 369 #  undef strtod
 370 # endif
 371 # if HAS_GRADUAL_UNDERFLOW_PROBLEM
 372 #  define SET_ERRNO_UPON_GRADUAL_UNDERFLOW(RESULT) \
 373     do                                                          \
 374       {                                                         \
 375         if ((RESULT) != 0 && (RESULT) < MIN && (RESULT) > -MIN) \
 376           errno = ERANGE;                                       \
 377       }                                                         \
 378     while (0)
 379 # else
 380 #  define SET_ERRNO_UPON_GRADUAL_UNDERFLOW(RESULT) (void)0
 381 # endif
 382 #else
 383 # undef STRTOD
 384 # define STRTOD(NPTR,ENDPTR) \
 385    parse_number (NPTR, 10, 10, 1, radixchar, 'e', ENDPTR)
 386 # define SET_ERRNO_UPON_GRADUAL_UNDERFLOW(RESULT) (void)0
 387 #endif
 388 /* From here on, STRTOD refers to the underlying implementation.  It needs
 389    to handle only finite unsigned decimal numbers with non-null ENDPTR.  */
 390 {
 391   char radixchar;
 392   bool negative = false;
 393
 394   /* The number so far.  */
 395   DOUBLE num;
 396
 397   const char *s = nptr;
 398   const char *end;
 399   char *endbuf;
 400   int saved_errno = errno;
 401
 402   radixchar = decimal_point_char ();
 403
 404   /* Eat whitespace.  */
 405   while (locale_isspace (*s))
 406     ++s;
 407
 408   /* Get the sign.  */
 409   negative = *s == '-';
 410   if (*s == '-' || *s == '+')
 411     ++s;
 412
 413   num = STRTOD (s, &endbuf);
 414   SET_ERRNO_UPON_GRADUAL_UNDERFLOW (num);
 415   end = endbuf;
 416
 417   if (c_isdigit (s[*s == radixchar]))
 418     {
 419       /* If a hex float was converted incorrectly, do it ourselves.
 420          If the string starts with "0x" but does not contain digits,
 421          consume the "0" ourselves.  If a hex float is followed by a
 422          'p' but no exponent, then adjust the end pointer.  */
 423       if (*s == '0' && c_tolower (s[1]) == 'x')
 424         {
 425           if (! c_isxdigit (s[2 + (s[2] == radixchar)]))
 426             {
 427               end = s + 1;
 428
 429               /* strtod() on z/OS returns ERANGE for "0x".  */
 430               errno = saved_errno;
 431             }
 432           else if (end <= s + 2)
 433             {
 434               num = parse_number (s + 2, 16, 2, 4, radixchar, 'p', &endbuf);
 435               end = endbuf;
 436             }
 437           else
 438             {
 439               const char *p = s + 2;
 440               while (p < end && c_tolower (*p) != 'p')
 441                 p++;
 442               if (p < end && ! c_isdigit (p[1 + (p[1] == '-' || p[1] == '+')]))
 443                 {
 444                   char *dup = strdup (s);
 445                   errno = saved_errno;
 446                   if (!dup)
 447                     {
 448                       /* Not really our day, is it.  Rounding errors are
 449                          better than outright failure.  */
 450                       num =
 451                         parse_number (s + 2, 16, 2, 4, radixchar, 'p', &endbuf);
 452                     }
 453                   else
 454                     {
 455                       dup[p - s] = '\0';
 456                       num = STRTOD (dup, &endbuf);
 457                       SET_ERRNO_UPON_GRADUAL_UNDERFLOW (num);
 458                       saved_errno = errno;
 459                       free (dup);
 460                       errno = saved_errno;
 461                     }
 462                   end = p;
 463                 }
 464             }
 465         }
 466       else
 467         {
 468           /* If "1e 1" was misparsed as 10.0 instead of 1.0, re-do the
 469              underlying STRTOD on a copy of the original string
 470              truncated to avoid the bug.  */
 471           const char *e = s + 1;
 472           while (e < end && c_tolower (*e) != 'e')
 473             e++;
 474           if (e < end && ! c_isdigit (e[1 + (e[1] == '-' || e[1] == '+')]))
 475             {
 476               char *dup = strdup (s);
 477               errno = saved_errno;
 478               if (!dup)
 479                 {
 480                   /* Not really our day, is it.  Rounding errors are
 481                      better than outright failure.  */
 482                   num = parse_number (s, 10, 10, 1, radixchar, 'e', &endbuf);
 483                 }
 484               else
 485                 {
 486                   dup[e - s] = '\0';
 487                   num = STRTOD (dup, &endbuf);
 488                   SET_ERRNO_UPON_GRADUAL_UNDERFLOW (num);
 489                   saved_errno = errno;
 490                   free (dup);
 491                   errno = saved_errno;
 492                 }
 493               end = e;
 494             }
 495         }
 496
 497       s = end;
 498     }
 499
 500   /* Check for infinities and NaNs.  */
 501   else if (c_tolower (*s) == 'i'
 502            && c_tolower (s[1]) == 'n'
 503            && c_tolower (s[2]) == 'f')
 504     {
 505       s += 3;
 506       if (c_tolower (*s) == 'i'
 507           && c_tolower (s[1]) == 'n'
 508           && c_tolower (s[2]) == 'i'
 509           && c_tolower (s[3]) == 't'
 510           && c_tolower (s[4]) == 'y')
 511         s += 5;
 512       num = HUGE_VAL;
 513       errno = saved_errno;
 514     }
 515   else if (c_tolower (*s) == 'n'
 516            && c_tolower (s[1]) == 'a'
 517            && c_tolower (s[2]) == 'n')
 518     {
 519       s += 3;
 520       if (*s == '(')
 521         {
 522           const char *p = s + 1;
 523           while (c_isalnum (*p))
 524             p++;
 525           if (*p == ')')
 526             s = p + 1;
 527         }
 528
 529       /* If the underlying implementation misparsed the NaN, assume
 530          its result is incorrect, and return a NaN.  Normally it's
 531          better to use the underlying implementation's result, since a
 532          nice implementation populates the bits of the NaN according
 533          to interpreting n-char-sequence as a hexadecimal number.  */
 534       if (s != end || num == num)
 535         num = NAN;
 536       errno = saved_errno;
 537     }
 538   else
 539     {
 540       /* No conversion could be performed.  */
 541       errno = EINVAL;
 542       s = nptr;
 543     }
 544
 545   if (endptr != NULL)
 546     *endptr = (char *) s;
 547   /* Special case -0.0, since at least ICC miscompiles negation.  We
 548      can't use copysign(), as that drags in -lm on some platforms.  */
 549   if (!num && negative)
 550     return minus_zero ();
 551   return negative ? -num : num;
 552 }