1 /* $NetBSD: strtod.c,v 1.18 2021/05/06 16:15:33 christos Exp $ */
2
3 /****************************************************************
4
5 The author of this software is David M. Gay.
6
7 Copyright (C) 1998-2001 by Lucent Technologies
8 All Rights Reserved
9
10 Permission to use, copy, modify, and distribute this software and
11 its documentation for any purpose and without fee is hereby
12 granted, provided that the above copyright notice appear in all
13 copies and that both that the copyright notice and this
14 permission notice and warranty disclaimer appear in supporting
15 documentation, and that the name of Lucent or any of its entities
16 not be used in advertising or publicity pertaining to
17 distribution of the software without specific, written prior
18 permission.
19
20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
27 THIS SOFTWARE.
28
29 ****************************************************************/
30
31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
32 * with " at " changed at "@" and " dot " changed to "."). */
33
34 #include "namespace.h"
35 #include "gdtoaimp.h"
36 #ifndef NO_FENV_H
37 #include <fenv.h>
38 #endif
39
40 #ifdef USE_LOCALE
41 #include <locale.h>
42 #include "setlocale_local.h"
43 #endif
44
45 #ifdef IEEE_Arith
46 #ifndef NO_IEEE_Scale
47 #define Avoid_Underflow
48 #undef tinytens
49 /* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */
50 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
51 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
52 9007199254740992.*9007199254740992.e-256
53 };
54 #endif
55 #endif
56
57 #ifdef Honor_FLT_ROUNDS
58 #undef Check_FLT_ROUNDS
59 #define Check_FLT_ROUNDS
60 #else
61 #define Rounding Flt_Rounds
62 #endif
63
64 #ifndef __HAVE_LONG_DOUBLE
65 __strong_alias(_strtold, strtod)
66 __weak_alias(strtold, _strtold)
67 __strong_alias(_strtold_l, strtod_l)
68 __weak_alias(strtold_l, _strtold_l)
69 #endif
70
71 #ifdef Avoid_Underflow /*{*/
72 static double
73 sulp
74 #ifdef KR_headers
75 (x, scale) U *x; int scale;
76 #else
77 (U *x, int scale)
78 #endif
79 {
80 U u;
81 double rv;
82 int i;
83
84 rv = ulp(x);
85 if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
86 return rv; /* Is there an example where i <= 0 ? */
87 word0(&u) = Exp_1 + (i << Exp_shift);
88 word1(&u) = 0;
89 return rv * u.d;
90 }
91 #endif /*}*/
92
93 #if __GNUC_PREREQ__(9, 3)
94 __attribute__((__optimize__("O0")))
95 #endif
96 static double
_int_strtod_l(CONST char * s00,char ** se,locale_t loc)97 _int_strtod_l(CONST char *s00, char **se, locale_t loc)
98 {
99 #ifdef Avoid_Underflow
100 int scale;
101 #endif
102 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
103 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
104 CONST char *s, *s0, *s1;
105 double aadj;
106 Long L;
107 U adj, aadj1, rv, rv0;
108 ULong y, z;
109 Bigint *bb = NULL, *bb1, *bd0;
110 Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */
111 #ifdef Avoid_Underflow
112 ULong Lsb, Lsb1;
113 #endif
114 #ifdef SET_INEXACT
115 int inexact, oldinexact;
116 #endif
117 #ifdef USE_LOCALE /*{{*/
118 char *decimalpoint = localeconv_l(loc)->decimal_point;
119 size_t dplen = strlen(decimalpoint);
120 #endif /*USE_LOCALE}}*/
121
122 #ifdef Honor_FLT_ROUNDS /*{*/
123 int Rounding;
124 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
125 Rounding = Flt_Rounds;
126 #else /*}{*/
127 Rounding = 1;
128 switch(fegetround()) {
129 case FE_TOWARDZERO: Rounding = 0; break;
130 case FE_UPWARD: Rounding = 2; break;
131 case FE_DOWNWARD: Rounding = 3;
132 }
133 #endif /*}}*/
134 #endif /*}*/
135
136 sign = nz0 = nz = decpt = 0;
137 dval(&rv) = 0.;
138 for(s = s00;;s++) switch(*s) {
139 case '-':
140 sign = 1;
141 /* FALLTHROUGH */
142 case '+':
143 if (*++s)
144 goto break2;
145 /* FALLTHROUGH */
146 case 0:
147 goto ret0;
148 case '\t':
149 case '\n':
150 case '\v':
151 case '\f':
152 case '\r':
153 case ' ':
154 continue;
155 default:
156 goto break2;
157 }
158 break2:
159 if (*s == '0') {
160 #ifndef NO_HEX_FP /*{*/
161 {
162 static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
163 Long expt;
164 ULong bits[2];
165 switch(s[1]) {
166 case 'x':
167 case 'X':
168 {
169 #ifdef Honor_FLT_ROUNDS
170 FPI fpi1 = fpi;
171 fpi1.rounding = Rounding;
172 #else
173 #define fpi1 fpi
174 #endif
175 switch((i = gethex(&s, &fpi1, &expt, &bb, sign, loc)) & STRTOG_Retmask) {
176 case STRTOG_NoNumber:
177 s = s00;
178 sign = 0;
179 /* FALLTHROUGH */
180 case STRTOG_Zero:
181 break;
182 default:
183 if (bb) {
184 copybits(bits, fpi.nbits, bb);
185 Bfree(bb);
186 }
187 ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i);
188 }}
189 goto ret;
190 }
191 }
192 #endif /*}*/
193 nz0 = 1;
194 while(*++s == '0') ;
195 if (!*s)
196 goto ret;
197 }
198 s0 = s;
199 y = z = 0;
200 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
201 if (nd < 9)
202 y = 10*y + c - '0';
203 else if (nd < DBL_DIG + 2)
204 z = 10*z + c - '0';
205 nd0 = nd;
206 #ifdef USE_LOCALE
207 if (c == *decimalpoint) {
208 for(i = 1; decimalpoint[i]; ++i)
209 if (s[i] != decimalpoint[i])
210 goto dig_done;
211 s += i;
212 c = *s;
213 #else
214 if (c == '.') {
215 c = *++s;
216 #endif
217 decpt = 1;
218 if (!nd) {
219 for(; c == '0'; c = *++s)
220 nz++;
221 if (c > '0' && c <= '9') {
222 s0 = s;
223 nf += nz;
224 nz = 0;
225 goto have_dig;
226 }
227 goto dig_done;
228 }
229 for(; c >= '0' && c <= '9'; c = *++s) {
230 have_dig:
231 nz++;
232 if (c -= '0') {
233 nf += nz;
234 for(i = 1; i < nz; i++)
235 if (nd++ < 9)
236 y *= 10;
237 else if (nd <= DBL_DIG + 2)
238 z *= 10;
239 if (nd++ < 9)
240 y = 10*y + c;
241 else if (nd <= DBL_DIG + 2)
242 z = 10*z + c;
243 nz = 0;
244 }
245 }
246 }/*}*/
247 dig_done:
248 e = 0;
249 if (c == 'e' || c == 'E') {
250 if (!nd && !nz && !nz0) {
251 goto ret0;
252 }
253 s00 = s;
254 esign = 0;
255 switch(c = *++s) {
256 case '-':
257 esign = 1;
258 /* FALLTHROUGH */
259 case '+':
260 c = *++s;
261 }
262 if (c >= '0' && c <= '9') {
263 while(c == '0')
264 c = *++s;
265 if (c > '0' && c <= '9') {
266 L = c - '0';
267 s1 = s;
268 while((c = *++s) >= '0' && c <= '9')
269 L = 10*L + c - '0';
270 if (s - s1 > 8 || L > 19999)
271 /* Avoid confusion from exponents
272 * so large that e might overflow.
273 */
274 e = 19999; /* safe for 16 bit ints */
275 else
276 e = (int)L;
277 if (esign)
278 e = -e;
279 }
280 else
281 e = 0;
282 }
283 else
284 s = s00;
285 }
286 if (!nd) {
287 if (!nz && !nz0) {
288 #ifdef INFNAN_CHECK
289 /* Check for Nan and Infinity */
290 ULong bits[2];
291 static CONST FPI fpinan = /* only 52 explicit bits */
292 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
293 if (!decpt)
294 switch(c) {
295 case 'i':
296 case 'I':
297 if (match(&s,"nf")) {
298 --s;
299 if (!match(&s,"inity"))
300 ++s;
301 word0(&rv) = 0x7ff00000;
302 word1(&rv) = 0;
303 goto ret;
304 }
305 break;
306 case 'n':
307 case 'N':
308 if (match(&s, "an")) {
309 #ifndef No_Hex_NaN
310 if (*s == '(' /*)*/
311 && hexnan(&s, &fpinan, bits)
312 == STRTOG_NaNbits) {
313 word0(&rv) = 0x7ff00000 | bits[1];
314 word1(&rv) = bits[0];
315 }
316 else {
317 #endif
318 word0(&rv) = NAN_WORD0;
319 word1(&rv) = NAN_WORD1;
320 #ifndef No_Hex_NaN
321 }
322 #endif
323 goto ret;
324 }
325 }
326 #endif /* INFNAN_CHECK */
327 ret0:
328 s = s00;
329 sign = 0;
330 }
331 goto ret;
332 }
333 e1 = e -= nf;
334
335 /* Now we have nd0 digits, starting at s0, followed by a
336 * decimal point, followed by nd-nd0 digits. The number we're
337 * after is the integer represented by those digits times
338 * 10**e */
339
340 if (!nd0)
341 nd0 = nd;
342 k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
343 dval(&rv) = y;
344 if (k > 9) {
345 #ifdef SET_INEXACT
346 if (k > DBL_DIG)
347 oldinexact = get_inexact();
348 #endif
349 dval(&rv) = tens[k - 9] * dval(&rv) + z;
350 }
351 bd0 = 0;
352 if (nd <= DBL_DIG
353 #ifndef RND_PRODQUOT
354 #ifndef Honor_FLT_ROUNDS
355 && Flt_Rounds == 1
356 #endif
357 #endif
358 ) {
359 if (!e)
360 goto ret;
361 #ifndef ROUND_BIASED_without_Round_Up
362 if (e > 0) {
363 if (e <= Ten_pmax) {
364 #ifdef VAX
365 goto vax_ovfl_check;
366 #else
367 #ifdef Honor_FLT_ROUNDS
368 /* round correctly FLT_ROUNDS = 2 or 3 */
369 if (sign) {
370 rv.d = -rv.d;
371 sign = 0;
372 }
373 #endif
374 /* rv = */ rounded_product(dval(&rv), tens[e]);
375 goto ret;
376 #endif
377 }
378 i = DBL_DIG - nd;
379 if (e <= Ten_pmax + i) {
380 /* A fancier test would sometimes let us do
381 * this for larger i values.
382 */
383 #ifdef Honor_FLT_ROUNDS
384 /* round correctly FLT_ROUNDS = 2 or 3 */
385 if (sign) {
386 rv.d = -rv.d;
387 sign = 0;
388 }
389 #endif
390 e -= i;
391 dval(&rv) *= tens[i];
392 #ifdef VAX
393 /* VAX exponent range is so narrow we must
394 * worry about overflow here...
395 */
396 vax_ovfl_check:
397 word0(&rv) -= P*Exp_msk1;
398 /* rv = */ rounded_product(dval(&rv), tens[e]);
399 if ((word0(&rv) & Exp_mask)
400 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
401 goto ovfl;
402 word0(&rv) += P*Exp_msk1;
403 #else
404 /* rv = */ rounded_product(dval(&rv), tens[e]);
405 #endif
406 goto ret;
407 }
408 }
409 #ifndef Inaccurate_Divide
410 else if (e >= -Ten_pmax) {
411 #ifdef Honor_FLT_ROUNDS
412 /* round correctly FLT_ROUNDS = 2 or 3 */
413 if (sign) {
414 rv.d = -rv.d;
415 sign = 0;
416 }
417 #endif
418 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
419 goto ret;
420 }
421 #endif
422 #endif /* ROUND_BIASED_without_Round_Up */
423 }
424 e1 += nd - k;
425
426 #ifdef IEEE_Arith
427 #ifdef SET_INEXACT
428 inexact = 1;
429 if (k <= DBL_DIG)
430 oldinexact = get_inexact();
431 #endif
432 #ifdef Avoid_Underflow
433 scale = 0;
434 #endif
435 #ifdef Honor_FLT_ROUNDS
436 if (Rounding >= 2) {
437 if (sign)
438 Rounding = Rounding == 2 ? 0 : 2;
439 else
440 if (Rounding != 2)
441 Rounding = 0;
442 }
443 #endif
444 #endif /*IEEE_Arith*/
445
446 /* Get starting approximation = rv * 10**e1 */
447
448 if (e1 > 0) {
449 if ( (i = e1 & 15) !=0)
450 dval(&rv) *= tens[i];
451 if (e1 &= ~15) {
452 if (e1 > DBL_MAX_10_EXP) {
453 ovfl:
454 /* Can't trust HUGE_VAL */
455 #ifdef IEEE_Arith
456 #ifdef Honor_FLT_ROUNDS
457 switch(Rounding) {
458 case 0: /* toward 0 */
459 case 3: /* toward -infinity */
460 word0(&rv) = Big0;
461 word1(&rv) = Big1;
462 break;
463 default:
464 word0(&rv) = Exp_mask;
465 word1(&rv) = 0;
466 }
467 #else /*Honor_FLT_ROUNDS*/
468 word0(&rv) = Exp_mask;
469 word1(&rv) = 0;
470 #endif /*Honor_FLT_ROUNDS*/
471 #ifdef SET_INEXACT
472 /* set overflow bit */
473 dval(&rv0) = 1e300;
474 dval(&rv0) *= dval(&rv0);
475 #endif
476 #else /*IEEE_Arith*/
477 word0(&rv) = Big0;
478 word1(&rv) = Big1;
479 #endif /*IEEE_Arith*/
480 range_err:
481 if (bd0) {
482 Bfree(bb);
483 Bfree(bd);
484 Bfree(bs);
485 Bfree(bd0);
486 Bfree(delta);
487 }
488 #ifndef NO_ERRNO
489 errno = ERANGE;
490 #endif
491 goto ret;
492 }
493 e1 = (unsigned int)e1 >> 4;
494 for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
495 if (e1 & 1)
496 dval(&rv) *= bigtens[j];
497 /* The last multiplication could overflow. */
498 word0(&rv) -= P*Exp_msk1;
499 dval(&rv) *= bigtens[j];
500 if ((z = word0(&rv) & Exp_mask)
501 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
502 goto ovfl;
503 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
504 /* set to largest number */
505 /* (Can't trust DBL_MAX) */
506 word0(&rv) = Big0;
507 word1(&rv) = Big1;
508 }
509 else
510 word0(&rv) += P*Exp_msk1;
511 }
512 }
513 else if (e1 < 0) {
514 e1 = -e1;
515 if ( (i = e1 & 15) !=0)
516 dval(&rv) /= tens[i];
517 if (e1 >>= 4) {
518 if (e1 >= 1 << n_bigtens)
519 goto undfl;
520 #ifdef Avoid_Underflow
521 if (e1 & Scale_Bit)
522 scale = 2*P;
523 for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
524 if (e1 & 1)
525 dval(&rv) *= tinytens[j];
526 if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
527 >> Exp_shift)) > 0) {
528 /* scaled rv is denormal; zap j low bits */
529 if (j >= 32) {
530 word1(&rv) = 0;
531 if (j >= 53)
532 word0(&rv) = (P+2)*Exp_msk1;
533 else
534 word0(&rv) &= 0xffffffffU << (j-32);
535 }
536 else
537 word1(&rv) &= 0xffffffffU << j;
538 }
539 #else
540 for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
541 if (e1 & 1)
542 dval(&rv) *= tinytens[j];
543 /* The last multiplication could underflow. */
544 dval(&rv0) = dval(&rv);
545 dval(&rv) *= tinytens[j];
546 if (!dval(&rv)) {
547 dval(&rv) = 2.*dval(&rv0);
548 dval(&rv) *= tinytens[j];
549 #endif
550 if (!dval(&rv)) {
551 undfl:
552 dval(&rv) = 0.;
553 #ifdef Honor_FLT_ROUNDS
554 if (Rounding == 2)
555 word1(&rv) = 1;
556 #endif
557 goto range_err;
558 }
559 #ifndef Avoid_Underflow
560 word0(&rv) = Tiny0;
561 word1(&rv) = Tiny1;
562 /* The refinement below will clean
563 * this approximation up.
564 */
565 }
566 #endif
567 }
568 }
569
570 /* Now the hard part -- adjusting rv to the correct value.*/
571
572 /* Put digits into bd: true value = bd * 10^e */
573
574 bd0 = s2b(s0, nd0, nd, y, dplen);
575 if (bd0 == NULL)
576 goto ovfl;
577
578 for(;;) {
579 bd = Balloc(bd0->k);
580 if (bd == NULL)
581 goto ovfl;
582 Bcopy(bd, bd0);
583 bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
584 if (bb == NULL)
585 goto ovfl;
586 bs = i2b(1);
587 if (bs == NULL)
588 goto ovfl;
589
590 if (e >= 0) {
591 bb2 = bb5 = 0;
592 bd2 = bd5 = e;
593 }
594 else {
595 bb2 = bb5 = -e;
596 bd2 = bd5 = 0;
597 }
598 if (bbe >= 0)
599 bb2 += bbe;
600 else
601 bd2 -= bbe;
602 bs2 = bb2;
603 #ifdef Honor_FLT_ROUNDS
604 if (Rounding != 1)
605 bs2++;
606 #endif
607 #ifdef Avoid_Underflow
608 Lsb = LSB;
609 Lsb1 = 0;
610 j = bbe - scale;
611 i = j + bbbits - 1; /* logb(rv) */
612 j = P + 1 - bbbits;
613 if (i < Emin) { /* denormal */
614 i = Emin - i;
615 j -= i;
616 if (i < 32)
617 Lsb <<= i;
618 else
619 Lsb1 = Lsb << (i-32);
620 }
621 #else /*Avoid_Underflow*/
622 #ifdef Sudden_Underflow
623 #ifdef IBM
624 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
625 #else
626 j = P + 1 - bbbits;
627 #endif
628 #else /*Sudden_Underflow*/
629 j = bbe;
630 i = j + bbbits - 1; /* logb(&rv) */
631 if (i < Emin) /* denormal */
632 j += P - Emin;
633 else
634 j = P + 1 - bbbits;
635 #endif /*Sudden_Underflow*/
636 #endif /*Avoid_Underflow*/
637 bb2 += j;
638 bd2 += j;
639 #ifdef Avoid_Underflow
640 bd2 += scale;
641 #endif
642 i = bb2 < bd2 ? bb2 : bd2;
643 if (i > bs2)
644 i = bs2;
645 if (i > 0) {
646 bb2 -= i;
647 bd2 -= i;
648 bs2 -= i;
649 }
650 if (bb5 > 0) {
651 bs = pow5mult(bs, bb5);
652 if (bs == NULL)
653 goto ovfl;
654 bb1 = mult(bs, bb);
655 if (bb1 == NULL)
656 goto ovfl;
657 Bfree(bb);
658 bb = bb1;
659 }
660 if (bb2 > 0) {
661 bb = lshift(bb, bb2);
662 if (bb == NULL)
663 goto ovfl;
664 }
665 if (bd5 > 0) {
666 bd = pow5mult(bd, bd5);
667 if (bd == NULL)
668 goto ovfl;
669 }
670 if (bd2 > 0) {
671 bd = lshift(bd, bd2);
672 if (bd == NULL)
673 goto ovfl;
674 }
675 if (bs2 > 0) {
676 bs = lshift(bs, bs2);
677 if (bs == NULL)
678 goto ovfl;
679 }
680 delta = diff(bb, bd);
681 if (delta == NULL)
682 goto ovfl;
683 dsign = delta->sign;
684 delta->sign = 0;
685 i = cmp(delta, bs);
686 #ifdef Honor_FLT_ROUNDS
687 if (Rounding != 1) {
688 if (i < 0) {
689 /* Error is less than an ulp */
690 if (!delta->x[0] && delta->wds <= 1) {
691 /* exact */
692 #ifdef SET_INEXACT
693 inexact = 0;
694 #endif
695 break;
696 }
697 if (Rounding) {
698 if (dsign) {
699 dval(&adj) = 1.;
700 goto apply_adj;
701 }
702 }
703 else if (!dsign) {
704 dval(&adj) = -1.;
705 if (!word1(&rv)
706 && !(word0(&rv) & Frac_mask)) {
707 y = word0(&rv) & Exp_mask;
708 #ifdef Avoid_Underflow
709 if (!scale || y > 2*P*Exp_msk1)
710 #else
711 if (y)
712 #endif
713 {
714 delta = lshift(delta,Log2P);
715 if (delta == NULL)
716 goto ovfl;
717 if (cmp(delta, bs) <= 0)
718 dval(&adj) = -0.5;
719 }
720 }
721 apply_adj:
722 #ifdef Avoid_Underflow
723 if (scale && (y = word0(&rv) & Exp_mask)
724 <= 2*P*Exp_msk1)
725 word0(&adj) += (2*P+1)*Exp_msk1 - y;
726 #else
727 #ifdef Sudden_Underflow
728 if ((word0(&rv) & Exp_mask) <=
729 P*Exp_msk1) {
730 word0(&rv) += P*Exp_msk1;
731 dval(&rv) += adj*ulp(&rv);
732 word0(&rv) -= P*Exp_msk1;
733 }
734 else
735 #endif /*Sudden_Underflow*/
736 #endif /*Avoid_Underflow*/
737 dval(&rv) += adj.d*ulp(&rv);
738 }
739 break;
740 }
741 dval(&adj) = ratio(delta, bs);
742 if (adj.d < 1.)
743 dval(&adj) = 1.;
744 if (adj.d <= 0x7ffffffe) {
745 /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
746 y = adj.d;
747 if (y != adj.d) {
748 if (!(((unsigned int)Rounding>>1) ^ (unsigned int)dsign))
749 y++;
750 dval(&adj) = y;
751 }
752 }
753 #ifdef Avoid_Underflow
754 if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
755 word0(&adj) += (2*P+1)*Exp_msk1 - y;
756 #else
757 #ifdef Sudden_Underflow
758 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
759 word0(&rv) += P*Exp_msk1;
760 dval(&adj) *= ulp(&rv);
761 if (dsign)
762 dval(&rv) += adj;
763 else
764 dval(&rv) -= adj;
765 word0(&rv) -= P*Exp_msk1;
766 goto cont;
767 }
768 #endif /*Sudden_Underflow*/
769 #endif /*Avoid_Underflow*/
770 dval(&adj) *= ulp(&rv);
771 if (dsign) {
772 if (word0(&rv) == Big0 && word1(&rv) == Big1)
773 goto ovfl;
774 dval(&rv) += adj.d;
775 }
776 else
777 dval(&rv) -= adj.d;
778 goto cont;
779 }
780 #endif /*Honor_FLT_ROUNDS*/
781
782 if (i < 0) {
783 /* Error is less than half an ulp -- check for
784 * special case of mantissa a power of two.
785 */
786 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
787 #ifdef IEEE_Arith
788 #ifdef Avoid_Underflow
789 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
790 #else
791 || (word0(&rv) & Exp_mask) <= Exp_msk1
792 #endif
793 #endif
794 ) {
795 #ifdef SET_INEXACT
796 if (!delta->x[0] && delta->wds <= 1)
797 inexact = 0;
798 #endif
799 break;
800 }
801 if (!delta->x[0] && delta->wds <= 1) {
802 /* exact result */
803 #ifdef SET_INEXACT
804 inexact = 0;
805 #endif
806 break;
807 }
808 delta = lshift(delta,Log2P);
809 if (delta == NULL)
810 goto ovfl;
811 if (cmp(delta, bs) > 0)
812 goto drop_down;
813 break;
814 }
815 if (i == 0) {
816 /* exactly half-way between */
817 if (dsign) {
818 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
819 && word1(&rv) == (
820 #ifdef Avoid_Underflow
821 (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
822 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
823 #endif
824 0xffffffff)) {
825 /*boundary case -- increment exponent*/
826 if (word0(&rv) == Big0 && word1(&rv) == Big1)
827 goto ovfl;
828 word0(&rv) = (word0(&rv) & Exp_mask)
829 + Exp_msk1
830 #ifdef IBM
831 | Exp_msk1 >> 4
832 #endif
833 ;
834 word1(&rv) = 0;
835 #ifdef Avoid_Underflow
836 dsign = 0;
837 #endif
838 break;
839 }
840 }
841 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
842 drop_down:
843 /* boundary case -- decrement exponent */
844 #ifdef Sudden_Underflow /*{{*/
845 L = word0(&rv) & Exp_mask;
846 #ifdef IBM
847 if (L < Exp_msk1)
848 #else
849 #ifdef Avoid_Underflow
850 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
851 #else
852 if (L <= Exp_msk1)
853 #endif /*Avoid_Underflow*/
854 #endif /*IBM*/
855 goto undfl;
856 L -= Exp_msk1;
857 #else /*Sudden_Underflow}{*/
858 #ifdef Avoid_Underflow
859 if (scale) {
860 L = word0(&rv) & Exp_mask;
861 if (L <= (2*P+1)*Exp_msk1) {
862 if (L > (P+2)*Exp_msk1)
863 /* round even ==> */
864 /* accept rv */
865 break;
866 /* rv = smallest denormal */
867 goto undfl;
868 }
869 }
870 #endif /*Avoid_Underflow*/
871 L = (word0(&rv) & Exp_mask) - Exp_msk1;
872 #endif /*Sudden_Underflow}}*/
873 word0(&rv) = L | Bndry_mask1;
874 word1(&rv) = 0xffffffff;
875 #ifdef IBM
876 goto cont;
877 #else
878 break;
879 #endif
880 }
881 #ifndef ROUND_BIASED
882 #ifdef Avoid_Underflow
883 if (Lsb1) {
884 if (!(word0(&rv) & Lsb1))
885 break;
886 }
887 else if (!(word1(&rv) & Lsb))
888 break;
889 #else
890 if (!(word1(&rv) & LSB))
891 break;
892 #endif
893 #endif
894 if (dsign)
895 #ifdef Avoid_Underflow
896 dval(&rv) += sulp(&rv, scale);
897 #else
898 dval(&rv) += ulp(&rv);
899 #endif
900 #ifndef ROUND_BIASED
901 else {
902 #ifdef Avoid_Underflow
903 dval(&rv) -= sulp(&rv, scale);
904 #else
905 dval(&rv) -= ulp(&rv);
906 #endif
907 #ifndef Sudden_Underflow
908 if (!dval(&rv))
909 goto undfl;
910 #endif
911 }
912 #ifdef Avoid_Underflow
913 dsign = 1 - dsign;
914 #endif
915 #endif
916 break;
917 }
918 if ((aadj = ratio(delta, bs)) <= 2.) {
919 if (dsign)
920 aadj = dval(&aadj1) = 1.;
921 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
922 #ifndef Sudden_Underflow
923 if (word1(&rv) == Tiny1 && !word0(&rv))
924 goto undfl;
925 #endif
926 aadj = 1.;
927 dval(&aadj1) = -1.;
928 }
929 else {
930 /* special case -- power of FLT_RADIX to be */
931 /* rounded down... */
932
933 if (aadj < 2./FLT_RADIX)
934 aadj = 1./FLT_RADIX;
935 else
936 aadj *= 0.5;
937 dval(&aadj1) = -aadj;
938 }
939 }
940 else {
941 aadj *= 0.5;
942 dval(&aadj1) = dsign ? aadj : -aadj;
943 #ifdef Check_FLT_ROUNDS
944 /* CONSTCOND */
945 switch(Rounding) {
946 case 2: /* towards +infinity */
947 dval(&aadj1) -= 0.5;
948 break;
949 case 0: /* towards 0 */
950 case 3: /* towards -infinity */
951 dval(&aadj1) += 0.5;
952 }
953 #else
954 /* CONSTCOND */
955 if (Flt_Rounds == 0)
956 dval(&aadj1) += 0.5;
957 #endif /*Check_FLT_ROUNDS*/
958 }
959 y = word0(&rv) & Exp_mask;
960
961 /* Check for overflow */
962
963 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
964 dval(&rv0) = dval(&rv);
965 word0(&rv) -= P*Exp_msk1;
966 dval(&adj) = dval(&aadj1) * ulp(&rv);
967 dval(&rv) += dval(&adj);
968 if ((word0(&rv) & Exp_mask) >=
969 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
970 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
971 goto ovfl;
972 word0(&rv) = Big0;
973 word1(&rv) = Big1;
974 goto cont;
975 }
976 else
977 word0(&rv) += P*Exp_msk1;
978 }
979 else {
980 #ifdef Avoid_Underflow
981 if (scale && y <= 2*P*Exp_msk1) {
982 if (aadj <= 0x7fffffff) {
983 if ((z = aadj) == 0)
984 z = 1;
985 aadj = z;
986 dval(&aadj1) = dsign ? aadj : -aadj;
987 }
988 word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
989 }
990 dval(&adj) = dval(&aadj1) * ulp(&rv);
991 dval(&rv) += dval(&adj);
992 #else
993 #ifdef Sudden_Underflow
994 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
995 dval(&rv0) = dval(&rv);
996 word0(&rv) += P*Exp_msk1;
997 dval(&adj) = dval(&aadj1) * ulp(&rv);
998 dval(&rv) += dval(&adj);
999 #ifdef IBM
1000 if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
1001 #else
1002 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
1003 #endif
1004 {
1005 if (word0(&rv0) == Tiny0
1006 && word1(&rv0) == Tiny1)
1007 goto undfl;
1008 word0(&rv) = Tiny0;
1009 word1(&rv) = Tiny1;
1010 goto cont;
1011 }
1012 else
1013 word0(&rv) -= P*Exp_msk1;
1014 }
1015 else {
1016 dval(&adj) = dval(&aadj1) * ulp(&rv);
1017 dval(&rv) += dval(&adj);
1018 }
1019 #else /*Sudden_Underflow*/
1020 /* Compute dval(&adj) so that the IEEE rounding rules will
1021 * correctly round rv + dval(&adj) in some half-way cases.
1022 * If rv * ulp(&rv) is denormalized (i.e.,
1023 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1024 * trouble from bits lost to denormalization;
1025 * example: 1.2e-307 .
1026 */
1027 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
1028 dval(&aadj1) = (double)(int)(aadj + 0.5);
1029 if (!dsign)
1030 dval(&aadj1) = -dval(&aadj1);
1031 }
1032 dval(&adj) = dval(&aadj1) * ulp(&rv);
1033 dval(&rv) += adj;
1034 #endif /*Sudden_Underflow*/
1035 #endif /*Avoid_Underflow*/
1036 }
1037 z = word0(&rv) & Exp_mask;
1038 #ifndef SET_INEXACT
1039 #ifdef Avoid_Underflow
1040 if (!scale)
1041 #endif
1042 if (y == z) {
1043 /* Can we stop now? */
1044 L = (Long)aadj;
1045 aadj -= L;
1046 /* The tolerances below are conservative. */
1047 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
1048 if (aadj < .4999999 || aadj > .5000001)
1049 break;
1050 }
1051 else if (aadj < .4999999/FLT_RADIX)
1052 break;
1053 }
1054 #endif
1055 cont:
1056 Bfree(bb);
1057 Bfree(bd);
1058 Bfree(bs);
1059 Bfree(delta);
1060 }
1061 Bfree(bb);
1062 Bfree(bd);
1063 Bfree(bs);
1064 Bfree(bd0);
1065 Bfree(delta);
1066 #ifdef SET_INEXACT
1067 if (inexact) {
1068 if (!oldinexact) {
1069 word0(&rv0) = Exp_1 + (70 << Exp_shift);
1070 word1(&rv0) = 0;
1071 dval(&rv0) += 1.;
1072 }
1073 }
1074 else if (!oldinexact)
1075 clear_inexact();
1076 #endif
1077 #ifdef Avoid_Underflow
1078 if (scale) {
1079 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
1080 word1(&rv0) = 0;
1081 dval(&rv) *= dval(&rv0);
1082 #ifndef NO_ERRNO
1083 /* try to avoid the bug of testing an 8087 register value */
1084 #ifdef IEEE_Arith
1085 if (!(word0(&rv) & Exp_mask))
1086 #else
1087 if (word0(&rv) == 0 && word1(&rv) == 0)
1088 #endif
1089 errno = ERANGE;
1090 #endif
1091 }
1092 #endif /* Avoid_Underflow */
1093 #ifdef SET_INEXACT
1094 if (inexact && !(word0(&rv) & Exp_mask)) {
1095 /* set underflow bit */
1096 dval(&rv0) = 1e-300;
1097 dval(&rv0) *= dval(&rv0);
1098 }
1099 #endif
1100 ret:
1101 if (se)
1102 *se = __UNCONST(s);
1103 return sign ? -dval(&rv) : dval(&rv);
1104 }
1105
1106 double
1107 strtod(CONST char *s, char **sp)
1108 {
1109 return _int_strtod_l(s, sp, _current_locale());
1110 }
1111
1112 #ifdef __weak_alias
1113 __weak_alias(strtod_l, _strtod_l)
1114 #endif
1115
1116 double
1117 strtod_l(CONST char *s, char **sp, locale_t loc)
1118 {
1119 return _int_strtod_l(s, sp, loc);
1120 }
1121