1 /*
2 FUNCTION
3 <<strtod>>, <<strtof>>, <<strtold>>, <<strtod_l>>, <<strtof_l>>, <<strtold_l>>---string to double or float
4
5 INDEX
6 strtod
7
8 INDEX
9 strtof
10
11 INDEX
12 strtold
13
14 INDEX
15 strtod_l
16
17 INDEX
18 strtof_l
19
20 INDEX
21 strtold_l
22
23 INDEX
24 _strtod_r
25
26 SYNOPSIS
27 #include <stdlib.h>
28 double strtod(const char *restrict <[str]>, char **restrict <[tail]>);
29 float strtof(const char *restrict <[str]>, char **restrict <[tail]>);
30 long double strtold(const char *restrict <[str]>,
31 char **restrict <[tail]>);
32
33 #include <stdlib.h>
34 double strtod_l(const char *restrict <[str]>, char **restrict <[tail]>,
35 locale_t <[locale]>);
36 float strtof_l(const char *restrict <[str]>, char **restrict <[tail]>,
37 locale_t <[locale]>);
38 long double strtold_l(const char *restrict <[str]>,
39 char **restrict <[tail]>,
40 locale_t <[locale]>);
41
42 double _strtod_r(void *<[reent]>,
43 const char *restrict <[str]>, char **restrict <[tail]>);
44
45 DESCRIPTION
46 <<strtod>>, <<strtof>>, <<strtold>> parse the character string
47 <[str]>, producing a substring which can be converted to a double,
48 float, or long double value, respectively. The substring converted
49 is the longest initial subsequence of <[str]>, beginning with the
50 first non-whitespace character, that has one of these formats:
51 .[+|-]<[digits]>[.[<[digits]>]][(e|E)[+|-]<[digits]>]
52 .[+|-].<[digits]>[(e|E)[+|-]<[digits]>]
53 .[+|-](i|I)(n|N)(f|F)[(i|I)(n|N)(i|I)(t|T)(y|Y)]
54 .[+|-](n|N)(a|A)(n|N)[<(>[<[hexdigits]>]<)>]
55 .[+|-]0(x|X)<[hexdigits]>[.[<[hexdigits]>]][(p|P)[+|-]<[digits]>]
56 .[+|-]0(x|X).<[hexdigits]>[(p|P)[+|-]<[digits]>]
57 The substring contains no characters if <[str]> is empty, consists
58 entirely of whitespace, or if the first non-whitespace
59 character is something other than <<+>>, <<->>, <<.>>, or a
60 digit, and cannot be parsed as infinity or NaN. If the platform
61 does not support NaN, then NaN is treated as an empty substring.
62 If the substring is empty, no conversion is done, and
63 the value of <[str]> is stored in <<*<[tail]>>>. Otherwise,
64 the substring is converted, and a pointer to the final string
65 (which will contain at least the terminating null character of
66 <[str]>) is stored in <<*<[tail]>>>. If you want no
67 assignment to <<*<[tail]>>>, pass a null pointer as <[tail]>.
68
69 This implementation returns the nearest machine number to the
70 input decimal string. Ties are broken by using the IEEE
71 round-even rule. However, <<strtof>> is currently subject to
72 double rounding errors.
73
74 <<strtod_l>>, <<strtof_l>>, <<strtold_l>> are like <<strtod>>,
75 <<strtof>>, <<strtold>> but perform the conversion based on the
76 locale specified by the locale object locale. If <[locale]> is
77 LC_GLOBAL_LOCALE or not a valid locale object, the behaviour is
78 undefined.
79
80 The alternate function <<_strtod_r>> is a reentrant version.
81 The extra argument <[reent]> is a pointer to a reentrancy structure.
82
83 RETURNS
84 These functions return the converted substring value, if any. If
85 no conversion could be performed, 0 is returned. If the correct
86 value is out of the range of representable values, plus or minus
87 <<HUGE_VAL>> (<<HUGE_VALF>>, <<HUGE_VALL>>) is returned, and
88 <<ERANGE>> is stored in errno. If the correct value would cause
89 underflow, 0 is returned and <<ERANGE>> is stored in errno.
90
91 PORTABILITY
92 <<strtod>> is ANSI.
93 <<strtof>>, <<strtold>> are C99.
94 <<strtod_l>>, <<strtof_l>>, <<strtold_l>> are GNU extensions.
95
96 Supporting OS subroutines required: <<close>>, <<fstat>>, <<isatty>>,
97 <<lseek>>, <<read>>, <<sbrk>>, <<write>>.
98 */
99
100 /****************************************************************
101
102 The author of this software is David M. Gay.
103
104 Copyright (C) 1998-2001 by Lucent Technologies
105 All Rights Reserved
106
107 Permission to use, copy, modify, and distribute this software and
108 its documentation for any purpose and without fee is hereby
109 granted, provided that the above copyright notice appear in all
110 copies and that both that the copyright notice and this
111 permission notice and warranty disclaimer appear in supporting
112 documentation, and that the name of Lucent or any of its entities
113 not be used in advertising or publicity pertaining to
114 distribution of the software without specific, written prior
115 permission.
116
117 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
118 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
119 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
120 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
121 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
122 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
123 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
124 THIS SOFTWARE.
125
126 ****************************************************************/
127
128 /* Please send bug reports to David M. Gay (dmg at acm dot org,
129 * with " at " changed at "@" and " dot " changed to "."). */
130
131 /* Original file gdtoa-strtod.c Modified 06-21-2006 by Jeff Johnston to work within newlib. */
132
133 #define _GNU_SOURCE
134 #include <errno.h>
135 #include <stdlib.h>
136 #include <string.h>
137 #include "mprec.h"
138 #include "gdtoa.h"
139 #include "../locale/setlocale.h"
140
141 /* #ifndef NO_FENV_H */
142 /* #include <fenv.h> */
143 /* #endif */
144
145 #include "locale.h"
146
147 #ifdef IEEE_Arith
148 #ifndef NO_IEEE_Scale
149 #define Avoid_Underflow
150 #undef tinytens
151 /* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */
152 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
153 static const double tinytens[] = { 1e-16, 1e-32,
154 #ifdef _DOUBLE_IS_32BITS
155 0.0, 0.0, 0.0
156 #else
157 1e-64, 1e-128,
158 9007199254740992. * 9007199254740992.e-256
159 #endif
160 };
161
162 #endif
163 #endif
164
165 #ifdef Honor_FLT_ROUNDS
166 #define Rounding rounding
167 #undef Check_FLT_ROUNDS
168 #define Check_FLT_ROUNDS
169 #else
170 #define Rounding Flt_Rounds
171 #endif
172
173 #ifdef IEEE_MC68k
174 #define _0 0
175 #define _1 1
176 #else
177 #define _0 1
178 #define _1 0
179 #endif
180
181 #ifdef Avoid_Underflow /*{*/
182 static double
sulp(U x,int scale)183 sulp (U x,
184 int scale)
185 {
186 U u;
187 double rv;
188 int i;
189
190 rv = ulp(dval(x));
191 if (!scale || (i = 2*P + 1 - ((dword0(x) & Exp_mask) >> Exp_shift)) <= 0)
192 return rv; /* Is there an example where i <= 0 ? */
193 dword0(u) = Exp_1 + ((__int32_t)i << Exp_shift);
194 #ifndef _DOUBLE_IS_32BITS
195 dword1(u) = 0;
196 #endif
197 return rv * u.d;
198 }
199 #endif /*}*/
200
201
202 #ifndef NO_HEX_FP
203
204 static void
ULtod(__ULong * L,__ULong * bits,Long exp,int k)205 ULtod (__ULong *L,
206 __ULong *bits,
207 Long exp,
208 int k)
209 {
210 switch(k & STRTOG_Retmask) {
211 case STRTOG_NoNumber:
212 case STRTOG_Zero:
213 L[0] = L[1] = 0;
214 break;
215
216 case STRTOG_Denormal:
217 L[_1] = bits[0];
218 L[_0] = bits[1];
219 break;
220
221 case STRTOG_Normal:
222 case STRTOG_NaNbits:
223 L[_1] = bits[0];
224 L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
225 break;
226
227 case STRTOG_Infinite:
228 L[_0] = 0x7ff00000;
229 L[_1] = 0;
230 break;
231
232 case STRTOG_NaN:
233 L[_0] = 0x7fffffff;
234 L[_1] = (__ULong)-1;
235 }
236 if (k & STRTOG_Neg)
237 L[_0] |= 0x80000000L;
238 }
239 #endif /* !NO_HEX_FP */
240
241 double
strtod_l(const char * __restrict s00,char ** __restrict se,locale_t loc)242 strtod_l (const char *__restrict s00, char **__restrict se,
243 locale_t loc)
244 {
245 #ifdef Avoid_Underflow
246 int scale;
247 #endif
248 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
249 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
250 const char *s, *s0, *s1;
251 double aadj, adj;
252 U aadj1, rv, rv0;
253 Long L;
254 __ULong y, z;
255 _Bigint *bb = NULL, *bb1, *bd = NULL, *bd0, *bs = NULL, *delta = NULL;
256 #ifdef Avoid_Underflow
257 __ULong Lsb, Lsb1;
258 #endif
259 #ifdef SET_INEXACT
260 int inexact, oldinexact;
261 #endif
262 #ifdef Honor_FLT_ROUNDS
263 int rounding;
264 #endif
265 #ifdef __HAVE_LOCALE_INFO__
266 const char *decimal_point = __get_numeric_locale(loc)->decimal_point;
267 const int dec_len = strlen (decimal_point);
268 #else
269 const char *decimal_point = ".";
270 const int dec_len = 1;
271
272 (void) loc;
273 #endif
274
275 delta = bs = bd = NULL;
276 sign = nz0 = nz = decpt = 0;
277 dval(rv) = 0.;
278 for(s = s00;;s++) switch(*s) {
279 case '-':
280 sign = 1;
281 __PICOLIBC_FALLTHROUGH;
282 case '+':
283 if (*++s)
284 goto break2;
285 __PICOLIBC_FALLTHROUGH;
286 case 0:
287 goto ret0;
288 case '\t':
289 case '\n':
290 case '\v':
291 case '\f':
292 case '\r':
293 case ' ':
294 continue;
295 default:
296 goto break2;
297 }
298 break2:
299 if (*s == '0') {
300 #ifndef NO_HEX_FP
301 {
302 static const FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
303 Long exp;
304 __ULong bits[2];
305 switch(s[1]) {
306 case 'x':
307 case 'X':
308 /* If the number is not hex, then the parse of
309 0 is still valid. */
310 s00 = s + 1;
311 {
312 #if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
313 FPI fpi1 = fpi;
314 switch(fegetround()) {
315 case FE_TOWARDZERO: fpi1.rounding = 0; break;
316 case FE_UPWARD: fpi1.rounding = 2; break;
317 case FE_DOWNWARD: fpi1.rounding = 3;
318 }
319 #else
320 #define fpi1 fpi
321 #endif
322 switch((i = gethex(&s, &fpi1, &exp, &bb, sign, loc)) & STRTOG_Retmask) {
323 case STRTOG_NoNumber:
324 s = s00;
325 sign = 0;
326 __PICOLIBC_FALLTHROUGH;
327 case STRTOG_Zero:
328 break;
329 default:
330 if (bb) {
331 copybits(bits, fpi.nbits, bb);
332 Bfree(bb);
333 }
334 ULtod(rv.i, bits, exp, i);
335 #ifndef NO_ERRNO
336 /* try to avoid the bug of testing an 8087 register value */
337 if ((dword0(rv)&Exp_mask) == 0)
338 errno = ERANGE;
339 #endif
340 }}
341 goto ret;
342 }
343 }
344 #endif
345 nz0 = 1;
346 while(*++s == '0') ;
347 if (!*s)
348 goto ret;
349 }
350 s0 = s;
351 y = z = 0;
352 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
353 if (nd < 9)
354 y = 10*y + c - '0';
355 else
356 z = 10*z + c - '0';
357 nd0 = nd;
358 if (strncmp (s, decimal_point, dec_len) == 0)
359 {
360 decpt = 1;
361 c = *(s += dec_len);
362 if (!nd) {
363 for(; c == '0'; c = *++s)
364 nz++;
365 if (c > '0' && c <= '9') {
366 s0 = s;
367 nf += nz;
368 nz = 0;
369 goto have_dig;
370 }
371 goto dig_done;
372 }
373 for(; c >= '0' && c <= '9'; c = *++s) {
374 have_dig:
375 nz++;
376 if (c -= '0') {
377 nf += nz;
378 for(i = 1; i < nz; i++)
379 if (nd++ < 9)
380 y *= 10;
381 else if (nd <= DBL_DIG + 1)
382 z *= 10;
383 if (nd++ < 9)
384 y = 10*y + c;
385 else if (nd <= DBL_DIG + 1)
386 z = 10*z + c;
387 nz = 0;
388 }
389 }
390 }
391 dig_done:
392 e = 0;
393 if (c == 'e' || c == 'E') {
394 if (!nd && !nz && !nz0) {
395 goto ret0;
396 }
397 s00 = s;
398 esign = 0;
399 switch(c = *++s) {
400 case '-':
401 esign = 1;
402 __PICOLIBC_FALLTHROUGH;
403 case '+':
404 c = *++s;
405 }
406 if (c >= '0' && c <= '9') {
407 while(c == '0')
408 c = *++s;
409 if (c > '0' && c <= '9') {
410 L = c - '0';
411 s1 = s;
412 while((c = *++s) >= '0' && c <= '9')
413 L = 10*L + c - '0';
414 if (s - s1 > 8 || L > 19999)
415 /* Avoid confusion from exponents
416 * so large that e might overflow.
417 */
418 e = 19999; /* safe for 16 bit ints */
419 else
420 e = (int)L;
421 if (esign)
422 e = -e;
423 }
424 else
425 e = 0;
426 }
427 else
428 s = s00;
429 }
430 if (!nd) {
431 if (!nz && !nz0) {
432 #ifdef INFNAN_CHECK
433 /* Check for Nan and Infinity */
434 __ULong bits[2];
435 static const FPI fpinan = /* only 52 explicit bits */
436 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
437 if (!decpt)
438 switch(c) {
439 case 'i':
440 case 'I':
441 if (match(&s,"nf")) {
442 --s;
443 if (!match(&s,"inity"))
444 ++s;
445 dword0(rv) = 0x7ff00000;
446 #ifndef _DOUBLE_IS_32BITS
447 dword1(rv) = 0;
448 #endif /*!_DOUBLE_IS_32BITS*/
449 goto ret;
450 }
451 break;
452 case 'n':
453 case 'N':
454 if (match(&s, "an")) {
455 #ifndef No_Hex_NaN
456 if (*s == '(' /*)*/
457 && hexnan(&s, &fpinan, bits)
458 == STRTOG_NaNbits) {
459 dword0(rv) = 0x7ff00000 | bits[1];
460 #ifndef _DOUBLE_IS_32BITS
461 dword1(rv) = bits[0];
462 #endif /*!_DOUBLE_IS_32BITS*/
463 }
464 else {
465 #endif
466 dval(rv) = nan ("");
467 #ifndef No_Hex_NaN
468 }
469 #endif
470 goto ret;
471 }
472 }
473 #endif /* INFNAN_CHECK */
474 ret0:
475 s = s00;
476 sign = 0;
477 }
478 goto ret;
479 }
480 e1 = e -= nf;
481
482 /* Now we have nd0 digits, starting at s0, followed by a
483 * decimal point, followed by nd-nd0 digits. The number we're
484 * after is the integer represented by those digits times
485 * 10**e */
486
487 if (!nd0)
488 nd0 = nd;
489 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
490 dval(rv) = y;
491 if (k > 9) {
492 #ifdef SET_INEXACT
493 if (k > DBL_DIG)
494 oldinexact = get_inexact();
495 #endif
496 dval(rv) = tens[k - 9] * dval(rv) + z;
497 }
498 bd0 = 0;
499 if (nd <= DBL_DIG
500 #ifndef RND_PRODQUOT
501 #ifndef Honor_FLT_ROUNDS
502 && Flt_Rounds == 1
503 #endif
504 #endif
505 ) {
506 if (!e)
507 goto ret;
508 if (e > 0) {
509 if (e <= Ten_pmax) {
510 #ifdef VAX
511 goto vax_ovfl_check;
512 #else
513 #ifdef Honor_FLT_ROUNDS
514 /* round correctly FLT_ROUNDS = 2 or 3 */
515 if (sign) {
516 dval(rv) = -dval(rv);
517 sign = 0;
518 }
519 #endif
520 /* rv = */ rounded_product(dval(rv), tens[e]);
521 goto ret;
522 #endif
523 }
524 i = DBL_DIG - nd;
525 if (e <= Ten_pmax + i) {
526 /* A fancier test would sometimes let us do
527 * this for larger i values.
528 */
529 #ifdef Honor_FLT_ROUNDS
530 /* round correctly FLT_ROUNDS = 2 or 3 */
531 if (sign) {
532 dval(rv) = -dval(rv);
533 sign = 0;
534 }
535 #endif
536 e -= i;
537 dval(rv) *= tens[i];
538 #ifdef VAX
539 /* VAX exponent range is so narrow we must
540 * worry about overflow here...
541 */
542 vax_ovfl_check:
543 dword0(rv) -= P*Exp_msk1;
544 /* rv = */ rounded_product(dval(rv), tens[e]);
545 if ((dword0(rv) & Exp_mask)
546 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
547 goto ovfl;
548 dword0(rv) += P*Exp_msk1;
549 #else
550 /* rv = */ rounded_product(dval(rv), tens[e]);
551 #endif
552 goto ret;
553 }
554 }
555 #ifndef Inaccurate_Divide
556 else if (e >= -Ten_pmax) {
557 #ifdef Honor_FLT_ROUNDS
558 /* round correctly FLT_ROUNDS = 2 or 3 */
559 if (sign) {
560 dval(rv) = -dval(rv);
561 sign = 0;
562 }
563 #endif
564 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
565 goto ret;
566 }
567 #endif
568 }
569 e1 += nd - k;
570
571 #ifdef IEEE_Arith
572 #ifdef SET_INEXACT
573 inexact = 1;
574 if (k <= DBL_DIG)
575 oldinexact = get_inexact();
576 #endif
577 #ifdef Avoid_Underflow
578 scale = 0;
579 #endif
580 #ifdef Honor_FLT_ROUNDS
581 if ((rounding = Flt_Rounds) >= 2) {
582 if (sign)
583 rounding = rounding == 2 ? 0 : 2;
584 else
585 if (rounding != 2)
586 rounding = 0;
587 }
588 #endif
589 #endif /*IEEE_Arith*/
590
591 /* Get starting approximation = rv * 10**e1 */
592
593 if (e1 > 0) {
594 if ( (i = e1 & 15) !=0)
595 dval(rv) *= tens[i];
596 if (e1 &= ~15) {
597 if (e1 > DBL_MAX_10_EXP) {
598 ovfl:
599 #ifndef NO_ERRNO
600 _REENT_ERRNO(ptr) = ERANGE;
601 #endif
602 /* Can't trust HUGE_VAL */
603 #ifdef IEEE_Arith
604 #ifdef Honor_FLT_ROUNDS
605 switch(rounding) {
606 case 0: /* toward 0 */
607 case 3: /* toward -infinity */
608 dword0(rv) = Big0;
609 #ifndef _DOUBLE_IS_32BITS
610 dword1(rv) = Big1;
611 #endif /*!_DOUBLE_IS_32BITS*/
612 break;
613 default:
614 dword0(rv) = Exp_mask;
615 #ifndef _DOUBLE_IS_32BITS
616 dword1(rv) = 0;
617 #endif /*!_DOUBLE_IS_32BITS*/
618 }
619 #else /*Honor_FLT_ROUNDS*/
620 dword0(rv) = Exp_mask;
621 #ifndef _DOUBLE_IS_32BITS
622 dword1(rv) = 0;
623 #endif /*!_DOUBLE_IS_32BITS*/
624 #endif /*Honor_FLT_ROUNDS*/
625 #ifdef SET_INEXACT
626 /* set overflow bit */
627 dval(rv0) = 1e300;
628 dval(rv0) *= dval(rv0);
629 #endif
630 #else /*IEEE_Arith*/
631 dword0(rv) = Big0;
632 #ifndef _DOUBLE_IS_32BITS
633 dword1(rv) = Big1;
634 #endif /*!_DOUBLE_IS_32BITS*/
635 #endif /*IEEE_Arith*/
636 if (bd0)
637 goto retfree;
638 goto ret;
639 }
640 e1 >>= 4;
641 for(j = 0; e1 > 1; j++, e1 >>= 1)
642 if (e1 & 1)
643 dval(rv) *= bigtens[j];
644 /* The last multiplication could overflow. */
645 dword0(rv) -= P*Exp_msk1;
646 dval(rv) *= bigtens[j];
647 if ((z = dword0(rv) & Exp_mask)
648 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
649 goto ovfl;
650 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
651 /* set to largest number */
652 /* (Can't trust DBL_MAX) */
653 dword0(rv) = Big0;
654 #ifndef _DOUBLE_IS_32BITS
655 dword1(rv) = Big1;
656 #endif /*!_DOUBLE_IS_32BITS*/
657 }
658 else
659 dword0(rv) += P*Exp_msk1;
660 }
661 }
662 else if (e1 < 0) {
663 e1 = -e1;
664 if ( (i = e1 & 15) !=0)
665 dval(rv) /= tens[i];
666 if (e1 >>= 4) {
667 if (e1 >= 1 << n_bigtens)
668 goto undfl;
669 #ifdef Avoid_Underflow
670 if (e1 & Scale_Bit)
671 scale = 2*P;
672 for(j = 0; e1 > 0; j++, e1 >>= 1)
673 if (e1 & 1)
674 dval(rv) *= tinytens[j];
675 if (scale && (j = 2*P + 1 - ((dword0(rv) & Exp_mask)
676 >> Exp_shift)) > 0) {
677 /* scaled rv is denormal; zap j low bits */
678 if (j >= 32) {
679 #ifndef _DOUBLE_IS_32BITS
680 dword1(rv) = 0;
681 #endif /*!_DOUBLE_IS_32BITS*/
682 if (j >= 53)
683 dword0(rv) = (P+2)*Exp_msk1;
684 else
685 dword0(rv) &= 0xffffffff << (j-32);
686 }
687 #ifndef _DOUBLE_IS_32BITS
688 else
689 dword1(rv) &= 0xffffffff << j;
690 #endif /*!_DOUBLE_IS_32BITS*/
691 }
692 #else
693 for(j = 0; e1 > 1; j++, e1 >>= 1)
694 if (e1 & 1)
695 dval(rv) *= tinytens[j];
696 /* The last multiplication could underflow. */
697 dval(rv0) = dval(rv);
698 dval(rv) *= tinytens[j];
699 if (!dval(rv)) {
700 dval(rv) = 2.*dval(rv0);
701 dval(rv) *= tinytens[j];
702 #endif
703 if (!dval(rv)) {
704 undfl:
705 dval(rv) = 0.;
706 #ifndef NO_ERRNO
707 _REENT_ERRNO(ptr) = ERANGE;
708 #endif
709 if (bd0)
710 goto retfree;
711 goto ret;
712 }
713 #ifndef Avoid_Underflow
714 #ifndef _DOUBLE_IS_32BITS
715 dword0(rv) = Tiny0;
716 dword1(rv) = Tiny1;
717 #else
718 dword0(rv) = Tiny1;
719 #endif /*_DOUBLE_IS_32BITS*/
720 /* The refinement below will clean
721 * this approximation up.
722 */
723 }
724 #endif
725 }
726 }
727
728 /* Now the hard part -- adjusting rv to the correct value.*/
729
730 /* Put digits into bd: true value = bd * 10^e */
731
732 bd0 = s2b(s0, nd0, nd, y);
733 if (bd0 == NULL)
734 goto ovfl;
735
736 for(;;) {
737 bd = Balloc(bd0->_k);
738 if (bd == NULL)
739 goto ovfl;
740 Bcopy(bd, bd0);
741 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
742 if (bb == NULL)
743 goto ovfl;
744 bs = i2b(1);
745 if (bs == NULL)
746 goto ovfl;
747
748 if (e >= 0) {
749 bb2 = bb5 = 0;
750 bd2 = bd5 = e;
751 }
752 else {
753 bb2 = bb5 = -e;
754 bd2 = bd5 = 0;
755 }
756 if (bbe >= 0)
757 bb2 += bbe;
758 else
759 bd2 -= bbe;
760 bs2 = bb2;
761 #ifdef Honor_FLT_ROUNDS
762 if (rounding != 1)
763 bs2++;
764 #endif
765 #ifdef Avoid_Underflow
766 Lsb = LSB;
767 Lsb1 = 0;
768 j = bbe - scale;
769 i = j + bbbits - 1; /* logb(rv) */
770 j = P + 1 - bbbits;
771 if (i < Emin) { /* denormal */
772 i = Emin - i;
773 j -= i;
774 if (i < 32)
775 Lsb <<= i;
776 else
777 Lsb1 = Lsb << (i-32);
778 }
779 #else /*Avoid_Underflow*/
780 #ifdef Sudden_Underflow
781 #ifdef IBM
782 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
783 #else
784 j = P + 1 - bbbits;
785 #endif
786 #else /*Sudden_Underflow*/
787 j = bbe;
788 i = j + bbbits - 1; /* logb(rv) */
789 if (i < Emin) /* denormal */
790 j += P - Emin;
791 else
792 j = P + 1 - bbbits;
793 #endif /*Sudden_Underflow*/
794 #endif /*Avoid_Underflow*/
795 bb2 += j;
796 bd2 += j;
797 #ifdef Avoid_Underflow
798 bd2 += scale;
799 #endif
800 i = bb2 < bd2 ? bb2 : bd2;
801 if (i > bs2)
802 i = bs2;
803 if (i > 0) {
804 bb2 -= i;
805 bd2 -= i;
806 bs2 -= i;
807 }
808 if (bb5 > 0) {
809 bs = pow5mult(bs, bb5);
810 if (bs == NULL)
811 goto ovfl;
812 bb1 = mult(bs, bb);
813 if (bb1 == NULL)
814 goto ovfl;
815 Bfree( bb);
816 bb = bb1;
817 }
818 if (bb2 > 0) {
819 bb = lshift(bb, bb2);
820 if (bb == NULL)
821 goto ovfl;
822 }
823 if (bd5 > 0) {
824 bd = pow5mult( bd, bd5);
825 if (bd == NULL)
826 goto ovfl;
827 }
828 if (bd2 > 0) {
829 bd = lshift(bd, bd2);
830 if (bd == NULL)
831 goto ovfl;
832 }
833 if (bs2 > 0) {
834 bs = lshift(bs, bs2);
835 if (bs == NULL)
836 goto ovfl;
837 }
838 delta = diff(bb, bd);
839 if (delta == NULL)
840 goto ovfl;
841 dsign = delta->_sign;
842 delta->_sign = 0;
843 i = cmp(delta, bs);
844 #ifdef Honor_FLT_ROUNDS
845 if (rounding != 1) {
846 if (i < 0) {
847 /* Error is less than an ulp */
848 if (!delta->_x[0] && delta->_wds <= 1) {
849 /* exact */
850 #ifdef SET_INEXACT
851 inexact = 0;
852 #endif
853 break;
854 }
855 if (rounding) {
856 if (dsign) {
857 adj = 1.;
858 goto apply_adj;
859 }
860 }
861 else if (!dsign) {
862 adj = -1.;
863 if (!dword1(rv)
864 && !(dword0(rv) & Frac_mask)) {
865 y = dword0(rv) & Exp_mask;
866 #ifdef Avoid_Underflow
867 if (!scale || y > 2*P*Exp_msk1)
868 #else
869 if (y)
870 #endif
871 {
872 delta = lshift(delta,Log2P);
873 if (cmp(delta, bs) <= 0)
874 adj = -0.5;
875 }
876 }
877 apply_adj:
878 #ifdef Avoid_Underflow
879 if (scale && (y = dword0(rv) & Exp_mask)
880 <= 2*P*Exp_msk1)
881 dword0(adj) += (2*P+1)*Exp_msk1 - y;
882 #else
883 #ifdef Sudden_Underflow
884 if ((dword0(rv) & Exp_mask) <=
885 P*Exp_msk1) {
886 dword0(rv) += P*Exp_msk1;
887 dval(rv) += adj*ulp(dval(rv));
888 dword0(rv) -= P*Exp_msk1;
889 }
890 else
891 #endif /*Sudden_Underflow*/
892 #endif /*Avoid_Underflow*/
893 dval(rv) += adj*ulp(dval(rv));
894 }
895 break;
896 }
897 adj = ratio(delta, bs);
898 if (adj < 1.)
899 adj = 1.;
900 if (adj <= 0x7ffffffe) {
901 /* adj = rounding ? ceil(adj) : floor(adj); */
902 y = adj;
903 if (y != adj) {
904 if (!((rounding>>1) ^ dsign))
905 y++;
906 adj = y;
907 }
908 }
909 #ifdef Avoid_Underflow
910 if (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
911 dword0(adj) += (2*P+1)*Exp_msk1 - y;
912 #else
913 #ifdef Sudden_Underflow
914 if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
915 dword0(rv) += P*Exp_msk1;
916 adj *= ulp(dval(rv));
917 if (dsign)
918 dval(rv) += adj;
919 else
920 dval(rv) -= adj;
921 dword0(rv) -= P*Exp_msk1;
922 goto cont;
923 }
924 #endif /*Sudden_Underflow*/
925 #endif /*Avoid_Underflow*/
926 adj *= ulp(dval(rv));
927 if (dsign) {
928 if (dword0(rv) == Big0 && dword1(rv) == Big1)
929 goto ovfl;
930 dval(rv) += adj;
931 else
932 dval(rv) -= adj;
933 goto cont;
934 }
935 #endif /*Honor_FLT_ROUNDS*/
936
937 if (i < 0) {
938 /* Error is less than half an ulp -- check for
939 * special case of mantissa a power of two.
940 */
941 if (dsign || dword1(rv) || dword0(rv) & Bndry_mask
942 #ifdef IEEE_Arith
943 #ifdef Avoid_Underflow
944 || (dword0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
945 #else
946 || (dword0(rv) & Exp_mask) <= Exp_msk1
947 #endif
948 #endif
949 ) {
950 #ifdef SET_INEXACT
951 if (!delta->x[0] && delta->wds <= 1)
952 inexact = 0;
953 #endif
954 break;
955 }
956 if (!delta->_x[0] && delta->_wds <= 1) {
957 /* exact result */
958 #ifdef SET_INEXACT
959 inexact = 0;
960 #endif
961 break;
962 }
963 delta = lshift(delta,Log2P);
964 if (cmp(delta, bs) > 0)
965 goto drop_down;
966 break;
967 }
968 if (i == 0) {
969 /* exactly half-way between */
970 if (dsign) {
971 if ((dword0(rv) & Bndry_mask1) == Bndry_mask1
972 && dword1(rv) == (
973 #ifdef Avoid_Underflow
974 (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
975 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
976 #endif
977 0xffffffff)) {
978 /*boundary case -- increment exponent*/
979 if (dword0(rv) == Big0 && dword1(rv) == Big1)
980 goto ovfl;
981 dword0(rv) = (dword0(rv) & Exp_mask)
982 + Exp_msk1
983 #ifdef IBM
984 | Exp_msk1 >> 4
985 #endif
986 ;
987 #ifndef _DOUBLE_IS_32BITS
988 dword1(rv) = 0;
989 #endif /*!_DOUBLE_IS_32BITS*/
990 #ifdef Avoid_Underflow
991 dsign = 0;
992 #endif
993 break;
994 }
995 }
996 else if (!(dword0(rv) & Bndry_mask) && !dword1(rv)) {
997 drop_down:
998 /* boundary case -- decrement exponent */
999 #ifdef Sudden_Underflow /*{{*/
1000 L = dword0(rv) & Exp_mask;
1001 #ifdef IBM
1002 if (L < Exp_msk1)
1003 #else
1004 #ifdef Avoid_Underflow
1005 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
1006 #else
1007 if (L <= Exp_msk1)
1008 #endif /*Avoid_Underflow*/
1009 #endif /*IBM*/
1010 goto undfl;
1011 L -= Exp_msk1;
1012 #else /*Sudden_Underflow}{*/
1013 #ifdef Avoid_Underflow
1014 if (scale) {
1015 L = dword0(rv) & Exp_mask;
1016 if (L <= (Long) ((2*P+1)*Exp_msk1)) {
1017 if (L > (Long)((P+2)*Exp_msk1))
1018 /* round even ==> */
1019 /* accept rv */
1020 break;
1021 /* rv = smallest denormal */
1022 goto undfl;
1023 }
1024 }
1025 #endif /*Avoid_Underflow*/
1026 L = (dword0(rv) & Exp_mask) - Exp_msk1;
1027 #endif /*Sudden_Underflow}*/
1028 dword0(rv) = L | Bndry_mask1;
1029 #ifndef _DOUBLE_IS_32BITS
1030 dword1(rv) = 0xffffffff;
1031 #endif /*!_DOUBLE_IS_32BITS*/
1032 #ifdef IBM
1033 goto cont;
1034 #else
1035 break;
1036 #endif
1037 }
1038 #ifndef ROUND_BIASED
1039 #ifdef Avoid_Underflow
1040 if (Lsb1) {
1041 if (!(dword0(rv) & Lsb1))
1042 break;
1043 }
1044 else if (!(dword1(rv) & Lsb))
1045 break;
1046 #else
1047 if (!(dword1(rv) & LSB))
1048 break;
1049 #endif
1050 #endif
1051 if (dsign)
1052 #ifdef Avoid_Underflow
1053 dval(rv) += sulp(rv, scale);
1054 #else
1055 dval(rv) += ulp(dval(rv));
1056 #endif
1057 #ifndef ROUND_BIASED
1058 else {
1059 #ifdef Avoid_Underflow
1060 dval(rv) -= sulp(rv, scale);
1061 #else
1062 dval(rv) -= ulp(dval(rv));
1063 #endif
1064 #ifndef Sudden_Underflow
1065 if (!dval(rv))
1066 goto undfl;
1067 #endif
1068 }
1069 #ifdef Avoid_Underflow
1070 dsign = 1 - dsign;
1071 #endif
1072 #endif
1073 break;
1074 }
1075 if ((aadj = ratio(delta, bs)) <= 2.) {
1076 if (dsign)
1077 aadj = dval(aadj1) = 1.;
1078 else if (dword1(rv) || dword0(rv) & Bndry_mask) {
1079 #ifndef Sudden_Underflow
1080 if (dword1(rv) == Tiny1 && !dword0(rv))
1081 goto undfl;
1082 #endif
1083 aadj = 1.;
1084 dval(aadj1) = -1.;
1085 }
1086 else {
1087 /* special case -- power of FLT_RADIX to be */
1088 /* rounded down... */
1089
1090 if (aadj < 2./FLT_RADIX)
1091 aadj = 1./FLT_RADIX;
1092 else
1093 aadj *= 0.5;
1094 dval(aadj1) = -aadj;
1095 }
1096 }
1097 else {
1098 aadj *= 0.5;
1099 dval(aadj1) = dsign ? aadj : -aadj;
1100 #ifdef Check_FLT_ROUNDS
1101 switch(Rounding) {
1102 case 2: /* towards +infinity */
1103 dval(aadj1) -= 0.5;
1104 break;
1105 case 0: /* towards 0 */
1106 case 3: /* towards -infinity */
1107 dval(aadj1) += 0.5;
1108 }
1109 #else
1110 if (Flt_Rounds == 0)
1111 dval(aadj1) += 0.5;
1112 #endif /*Check_FLT_ROUNDS*/
1113 }
1114 y = dword0(rv) & Exp_mask;
1115
1116 /* Check for overflow */
1117
1118 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1119 dval(rv0) = dval(rv);
1120 dword0(rv) -= P*Exp_msk1;
1121 adj = dval(aadj1) * ulp(dval(rv));
1122 dval(rv) += adj;
1123 if ((dword0(rv) & Exp_mask) >=
1124 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1125 if (dword0(rv0) == Big0 && dword1(rv0) == Big1)
1126 goto ovfl;
1127 dword0(rv) = Big0;
1128 #ifndef _DOUBLE_IS_32BITS
1129 dword1(rv) = Big1;
1130 #endif /*!_DOUBLE_IS_32BITS*/
1131 goto cont;
1132 }
1133 else
1134 dword0(rv) += P*Exp_msk1;
1135 }
1136 else {
1137 #ifdef Avoid_Underflow
1138 if (scale && y <= 2*P*Exp_msk1) {
1139 if (aadj <= 0x7fffffff) {
1140 if ((z = aadj) == 0)
1141 z = 1;
1142 aadj = z;
1143 dval(aadj1) = dsign ? aadj : -aadj;
1144 }
1145 dword0(aadj1) += (2*P+1)*Exp_msk1 - y;
1146 }
1147 adj = dval(aadj1) * ulp(dval(rv));
1148 dval(rv) += adj;
1149 #else
1150 #ifdef Sudden_Underflow
1151 if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
1152 dval(rv0) = dval(rv);
1153 dword0(rv) += P*Exp_msk1;
1154 adj = dval(aadj1) * ulp(dval(rv));
1155 dval(rv) += adj;
1156 #ifdef IBM
1157 if ((dword0(rv) & Exp_mask) < P*Exp_msk1)
1158 #else
1159 if ((dword0(rv) & Exp_mask) <= P*Exp_msk1)
1160 #endif
1161 {
1162 if (dword0(rv0) == Tiny0
1163 && dword1(rv0) == Tiny1)
1164 goto undfl;
1165 #ifndef _DOUBLE_IS_32BITS
1166 dword0(rv) = Tiny0;
1167 dword1(rv) = Tiny1;
1168 #else
1169 dword0(rv) = Tiny1;
1170 #endif /*_DOUBLE_IS_32BITS*/
1171 goto cont;
1172 }
1173 else
1174 dword0(rv) -= P*Exp_msk1;
1175 }
1176 else {
1177 adj = dval(aadj1) * ulp(dval(rv));
1178 dval(rv) += adj;
1179 }
1180 #else /*Sudden_Underflow*/
1181 /* Compute adj so that the IEEE rounding rules will
1182 * correctly round rv + adj in some half-way cases.
1183 * If rv * ulp(rv) is denormalized (i.e.,
1184 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1185 * trouble from bits lost to denormalization;
1186 * example: 1.2e-307 .
1187 */
1188 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
1189 dval(aadj1) = (double)(int)(aadj + 0.5);
1190 if (!dsign)
1191 dval(aadj1) = -dval(aadj1);
1192 }
1193 adj = dval(aadj1) * ulp(dval(rv));
1194 dval(rv) += adj;
1195 #endif /*Sudden_Underflow*/
1196 #endif /*Avoid_Underflow*/
1197 }
1198 z = dword0(rv) & Exp_mask;
1199 #ifndef SET_INEXACT
1200 #ifdef Avoid_Underflow
1201 if (!scale)
1202 #endif
1203 if (y == z) {
1204 /* Can we stop now? */
1205 #ifndef _DOUBLE_IS_32BITS
1206 /* If FE_INVALID floating point exceptions are
1207 enabled, a conversion to a 32 bit value is
1208 dangerous. A positive double value can result
1209 in a negative 32 bit int, thus raising SIGFPE.
1210 To avoid this, always convert into 64 bit here. */
1211 __int64_t L = (__int64_t)aadj;
1212 #else
1213 L = (Long)aadj;
1214 #endif
1215 aadj -= L;
1216 /* The tolerances below are conservative. */
1217 if (dsign || dword1(rv) || dword0(rv) & Bndry_mask) {
1218 if (aadj < .4999999 || aadj > .5000001)
1219 break;
1220 }
1221 else if (aadj < .4999999/FLT_RADIX)
1222 break;
1223 }
1224 #endif
1225 cont:
1226 Bfree(bb);
1227 Bfree(bd);
1228 Bfree(bs);
1229 Bfree(delta);
1230 }
1231 #ifdef SET_INEXACT
1232 if (inexact) {
1233 if (!oldinexact) {
1234 dword0(rv0) = Exp_1 + (70 << Exp_shift);
1235 #ifndef _DOUBLE_IS_32BITS
1236 dword1(rv0) = 0;
1237 #endif /*!_DOUBLE_IS_32BITS*/
1238 dval(rv0) += 1.;
1239 }
1240 }
1241 else if (!oldinexact)
1242 clear_inexact();
1243 #endif
1244 #ifdef Avoid_Underflow
1245 if (scale) {
1246 dword0(rv0) = Exp_1 - 2*P*Exp_msk1;
1247 #ifndef _DOUBLE_IS_32BITS
1248 dword1(rv0) = 0;
1249 #endif /*!_DOUBLE_IS_32BITS*/
1250 dval(rv) *= dval(rv0);
1251 #ifndef NO_ERRNO
1252 /* try to avoid the bug of testing an 8087 register value */
1253 if ((dword0(rv) & Exp_mask) == 0)
1254 _REENT_ERRNO(ptr) = ERANGE;
1255 #endif
1256 }
1257 #endif /* Avoid_Underflow */
1258 #ifdef SET_INEXACT
1259 if (inexact && !(dword0(rv) & Exp_mask)) {
1260 /* set underflow bit */
1261 dval(rv0) = 1e-300;
1262 dval(rv0) *= dval(rv0);
1263 }
1264 #endif
1265 retfree:
1266 Bfree(bb);
1267 Bfree(bd);
1268 Bfree(bs);
1269 Bfree(bd0);
1270 Bfree(delta);
1271 ret:
1272 if (se)
1273 *se = (char *)s;
1274 return sign ? -dval(rv) : dval(rv);
1275 }
1276
1277 double
1278 strtod (const char *__restrict s00,
1279 char **__restrict se)
1280 {
1281 return strtod_l (s00, se, __get_current_locale ());
1282 }
1283
1284 #if defined(_HAVE_LONG_DOUBLE) && defined(_LDBL_EQ_DBL)
1285 #ifdef _HAVE_ALIAS_ATTRIBUTE
1286 #pragma GCC diagnostic ignored "-Wpragmas"
1287 #pragma GCC diagnostic ignored "-Wunknown-warning-option"
1288 #pragma GCC diagnostic ignored "-Wattribute-alias="
1289 extern long double strtold(const char *, char **) __attribute__ ((__alias__ ("strtod")));
1290 #else
1291 long double
1292 strtold (const char * nptr, char ** endptr)
1293 {
1294 return (long double) strtod(nptr, endptr);
1295 }
1296 #endif
1297 #endif
1298
1299 /*
1300 * These two functions are not quite correct as they return true for
1301 * zero, however they are 'good enough' for the test in strtof below
1302 * as we only need to know whether the double test is false when
1303 * the float test is true.
1304 */
1305 static inline int
1306 isdenorm(double d)
1307 {
1308 U u;
1309 dval(u) = d;
1310 return (dword0(u) & Exp_mask) == 0;
1311 }
1312
1313 static inline int
1314 isdenormf(float f)
1315 {
1316 union { float f; __uint32_t i; } u;
1317 u.f = f;
1318 return (u.i & 0x7f800000) == 0;
1319 }
1320
1321 float
1322 strtof_l (const char *__restrict s00, char **__restrict se, locale_t loc)
1323 {
1324 double val = strtod_l (s00, se, loc);
1325 if (isnan (val))
1326 return signbit (val) ? -nanf ("") : nanf ("");
1327 float retval = (float) val;
1328 #ifndef NO_ERRNO
1329 if ((isinf (retval) && !isinf (val)) || (isdenormf(retval) && !isdenorm(val)))
1330 _REENT_ERRNO(_REENT) = ERANGE;
1331 #endif
1332 return retval;
1333 }
1334
1335 float
1336 strtof (const char *__restrict s00,
1337 char **__restrict se)
1338 {
1339 double val = strtod_l (s00, se, __get_current_locale ());
1340 if (isnan (val))
1341 return signbit (val) ? -nanf ("") : nanf ("");
1342 float retval = (float) val;
1343 #ifndef NO_ERRNO
1344 if ((isinf (retval) && !isinf (val)) || (isdenormf(retval) && !isdenorm(val)))
1345 _REENT_ERRNO(_REENT) = ERANGE;
1346 #endif
1347 return retval;
1348 }
1349