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                         FALLTHROUGH;
283 		case '+':
284 			if (*++s)
285 				goto break2;
286                         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 				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                                 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