strtod: work around IRIX 6.5 bug
[gnulib.git] / lib / strtod.c
1 /* Copyright (C) 1991-1992, 1997, 1999, 2003, 2006, 2008-2010 Free Software
2    Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation; either version 3 of the License, or
7    (at your option) any later version.
8
9    This program 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 General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
16
17 #include <config.h>
18
19 #include <stdlib.h>
20
21 #include <ctype.h>
22 #include <errno.h>
23 #include <float.h>
24 #include <limits.h>
25 #include <math.h>
26 #include <stdbool.h>
27 #include <string.h>
28
29 #include "c-ctype.h"
30
31 #ifndef HAVE_LDEXP_IN_LIBC
32 #define HAVE_LDEXP_IN_LIBC 0
33 #endif
34 #ifndef HAVE_RAW_DECL_STRTOD
35 #define HAVE_RAW_DECL_STRTOD 0
36 #endif
37
38 /* Return true if C is a space in the current locale, avoiding
39    problems with signed char and isspace.  */
40 static bool
41 locale_isspace (char c)
42 {
43   unsigned char uc = c;
44   return isspace (uc) != 0;
45 }
46
47 #if !HAVE_LDEXP_IN_LIBC
48  #define ldexp dummy_ldexp
49  /* A dummy definition that will never be invoked.  */
50  static double ldexp (double x _GL_UNUSED, int exponent _GL_UNUSED)
51  {
52    abort ();
53    return 0.0;
54  }
55 #endif
56
57 /* Return X * BASE**EXPONENT.  Return an extreme value and set errno
58    to ERANGE if underflow or overflow occurs.  */
59 static double
60 scale_radix_exp (double x, int radix, long int exponent)
61 {
62   /* If RADIX == 10, this code is neither precise nor fast; it is
63      merely a straightforward and relatively portable approximation.
64      If N == 2, this code is precise on a radix-2 implementation,
65      albeit perhaps not fast if ldexp is not in libc.  */
66
67   long int e = exponent;
68
69   if (HAVE_LDEXP_IN_LIBC && radix == 2)
70     return ldexp (x, e < INT_MIN ? INT_MIN : INT_MAX < e ? INT_MAX : e);
71   else
72     {
73       double r = x;
74
75       if (r != 0)
76         {
77           if (e < 0)
78             {
79               while (e++ != 0)
80                 {
81                   r /= radix;
82                   if (r == 0 && x != 0)
83                     {
84                       errno = ERANGE;
85                       break;
86                     }
87                 }
88             }
89           else
90             {
91               while (e-- != 0)
92                 {
93                   if (r < -DBL_MAX / radix)
94                     {
95                       errno = ERANGE;
96                       return -HUGE_VAL;
97                     }
98                   else if (DBL_MAX / radix < r)
99                     {
100                       errno = ERANGE;
101                       return HUGE_VAL;
102                     }
103                   else
104                     r *= radix;
105                 }
106             }
107         }
108
109       return r;
110     }
111 }
112
113 /* Parse a number at NPTR; this is a bit like strtol (NPTR, ENDPTR)
114    except there are no leading spaces or signs or "0x", and ENDPTR is
115    nonnull.  The number uses a base BASE (either 10 or 16) fraction, a
116    radix RADIX (either 10 or 2) exponent, and exponent character
117    EXPCHAR.  To convert from a number of digits to a radix exponent,
118    multiply by RADIX_MULTIPLIER (either 1 or 4).  */
119 static double
120 parse_number (const char *nptr,
121               int base, int radix, int radix_multiplier, char expchar,
122               char **endptr)
123 {
124   const char *s = nptr;
125   bool got_dot = false;
126   long int exponent = 0;
127   double num = 0;
128
129   for (;; ++s)
130     {
131       int digit;
132       if (c_isdigit (*s))
133         digit = *s - '0';
134       else if (base == 16 && c_isxdigit (*s))
135         digit = c_tolower (*s) - ('a' - 10);
136       else if (! got_dot && *s == '.')
137         {
138           /* Record that we have found the decimal point.  */
139           got_dot = true;
140           continue;
141         }
142       else
143         /* Any other character terminates the number.  */
144         break;
145
146       /* Make sure that multiplication by base will not overflow.  */
147       if (num <= DBL_MAX / base)
148         num = num * base + digit;
149       else
150         {
151           /* The value of the digit doesn't matter, since we have already
152              gotten as many digits as can be represented in a `double'.
153              This doesn't necessarily mean the result will overflow.
154              The exponent may reduce it to within range.
155
156              We just need to record that there was another
157              digit so that we can multiply by 10 later.  */
158           exponent += radix_multiplier;
159         }
160
161       /* Keep track of the number of digits after the decimal point.
162          If we just divided by base here, we might lose precision.  */
163       if (got_dot)
164         exponent -= radix_multiplier;
165     }
166
167   if (c_tolower (*s) == expchar && ! locale_isspace (s[1]))
168     {
169       /* Add any given exponent to the implicit one.  */
170       int save = errno;
171       char *end;
172       long int value = strtol (s + 1, &end, 10);
173       errno = save;
174
175       if (s + 1 != end)
176         {
177           /* Skip past the exponent, and add in the implicit exponent,
178              resulting in an extreme value on overflow.  */
179           s = end;
180           exponent =
181             (exponent < 0
182              ? (value < LONG_MIN - exponent ? LONG_MIN : exponent + value)
183              : (LONG_MAX - exponent < value ? LONG_MAX : exponent + value));
184         }
185     }
186
187   *endptr = (char *) s;
188   return scale_radix_exp (num, radix, exponent);
189 }
190
191 static double underlying_strtod (const char *, char **);
192
193 /* Convert NPTR to a double.  If ENDPTR is not NULL, a pointer to the
194    character after the last one used in the number is put in *ENDPTR.  */
195 double
196 strtod (const char *nptr, char **endptr)
197 {
198   bool negative = false;
199
200   /* The number so far.  */
201   double num;
202
203   const char *s = nptr;
204   const char *end;
205   char *endbuf;
206   int saved_errno;
207
208   /* Eat whitespace.  */
209   while (locale_isspace (*s))
210     ++s;
211
212   /* Get the sign.  */
213   negative = *s == '-';
214   if (*s == '-' || *s == '+')
215     ++s;
216
217   saved_errno = errno;
218   num = underlying_strtod (s, &endbuf);
219   end = endbuf;
220
221   if (c_isdigit (s[*s == '.']))
222     {
223       /* If a hex float was converted incorrectly, do it ourselves.
224          If the string starts with "0x" but does not contain digits,
225          consume the "0" ourselves.  If a hex float is followed by a
226          'p' but no exponent, then adjust the end pointer.  */
227       if (*s == '0' && c_tolower (s[1]) == 'x')
228         {
229           if (! c_isxdigit (s[2 + (s[2] == '.')]))
230             end = s + 1;
231           else if (end <= s + 2)
232             {
233               num = parse_number (s + 2, 16, 2, 4, 'p', &endbuf);
234               end = endbuf;
235             }
236           else
237             {
238               const char *p = s + 2;
239               while (p < end && c_tolower (*p) != 'p')
240                 p++;
241               if (p < end && ! c_isdigit (p[1 + (p[1] == '-' || p[1] == '+')]))
242                 end = p;
243             }
244         }
245       else
246         {
247           /* If "1e 1" was misparsed as 10.0 instead of 1.0, re-do the
248              underlying strtod on a copy of the original string
249              truncated to avoid the bug.  */
250           const char *e = s + 1;
251           while (e < end && c_tolower (*e) != 'e')
252             e++;
253           if (e < end && ! c_isdigit (e[1 + (e[1] == '-' || e[1] == '+')]))
254             {
255               char *dup = strdup (s);
256               errno = saved_errno;
257               if (!dup)
258                 {
259                   /* Not really our day, is it.  Rounding errors are
260                      better than outright failure.  */
261                   num = parse_number (s, 10, 10, 1, 'e', &endbuf);
262                 }
263               else
264                 {
265                   dup[e - s] = '\0';
266                   num = underlying_strtod (dup, &endbuf);
267                   saved_errno = errno;
268                   free (dup);
269                   errno = saved_errno;
270                 }
271               end = e;
272             }
273         }
274
275       s = end;
276     }
277
278   /* Check for infinities and NaNs.  */
279   else if (c_tolower (*s) == 'i'
280            && c_tolower (s[1]) == 'n'
281            && c_tolower (s[2]) == 'f')
282     {
283       s += 3;
284       if (c_tolower (*s) == 'i'
285           && c_tolower (s[1]) == 'n'
286           && c_tolower (s[2]) == 'i'
287           && c_tolower (s[3]) == 't'
288           && c_tolower (s[4]) == 'y')
289         s += 5;
290       num = HUGE_VAL;
291     }
292   else if (c_tolower (*s) == 'n'
293            && c_tolower (s[1]) == 'a'
294            && c_tolower (s[2]) == 'n')
295     {
296       s += 3;
297       if (*s == '(')
298         {
299           const char *p = s + 1;
300           while (c_isalnum (*p))
301             p++;
302           if (*p == ')')
303             s = p + 1;
304         }
305
306       /* If the underlying implementation misparsed the NaN, assume
307          its result is incorrect, and return a NaN.  Normally it's
308          better to use the underlying implementation's result, since a
309          nice implementation populates the bits of the NaN according
310          to interpreting n-char-sequence as a hexadecimal number.  */
311       if (s != end)
312         num = NAN;
313     }
314   else
315     {
316       /* No conversion could be performed.  */
317       errno = EINVAL;
318       s = nptr;
319     }
320
321   if (endptr != NULL)
322     *endptr = (char *) s;
323   return negative ? -num : num;
324 }
325
326 /* The "underlying" strtod implementation.  This must be defined
327    after strtod because it #undefs strtod.  */
328 static double
329 underlying_strtod (const char *nptr, char **endptr)
330 {
331   if (HAVE_RAW_DECL_STRTOD)
332     {
333       /* Prefer the native strtod if available.  Usually it should
334          work and it should give more-accurate results than our
335          approximation.  */
336       #undef strtod
337       return strtod (nptr, endptr);
338     }
339   else
340     {
341       /* Approximate strtod well enough for this module.  There's no
342          need to handle anything but finite unsigned decimal
343          numbers with nonnull ENDPTR.  */
344       return parse_number (nptr, 10, 10, 1, 'e', endptr);
345     }
346 }