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