1 /* $NetBSD: strtod.c,v 1.20 2024/01/23 15:31:58 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 static double
_int_strtod_l(CONST char * s00,char ** se,locale_t loc)94 _int_strtod_l(CONST char *s00, char **se, locale_t loc)
95 {
96 #ifdef Avoid_Underflow
97 int scale;
98 #endif
99 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
100 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
101 CONST char *s, *s0, *s1;
102 double aadj;
103 Long L;
104 U adj, aadj1, rv, rv0;
105 ULong y, z;
106 Bigint *bb = NULL, *bb1, *bd0;
107 Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */
108 #ifdef Avoid_Underflow
109 ULong Lsb, Lsb1;
110 #endif
111 #ifdef SET_INEXACT
112 int inexact, oldinexact;
113 #endif
114 #ifdef USE_LOCALE /*{{*/
115 char *decimalpoint = localeconv_l(loc)->decimal_point;
116 size_t dplen = strlen(decimalpoint);
117 #endif /*USE_LOCALE}}*/
118
119 #ifdef Honor_FLT_ROUNDS /*{*/
120 int Rounding, decpt = 0;
121 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
122 Rounding = Flt_Rounds;
123 #else /*}{*/
124 Rounding = 1;
125 switch(fegetround()) {
126 case FE_TOWARDZERO: Rounding = 0; break;
127 case FE_UPWARD: Rounding = 2; break;
128 case FE_DOWNWARD: Rounding = 3;
129 }
130 #endif /*}}*/
131 #endif /*}*/
132
133 sign = nz0 = nz = 0;
134 dval(&rv) = 0.;
135 for(s = s00;;s++) switch(*s) {
136 case '-':
137 sign = 1;
138 /* FALLTHROUGH */
139 case '+':
140 if (*++s)
141 goto break2;
142 /* FALLTHROUGH */
143 case 0:
144 goto ret0;
145 case '\t':
146 case '\n':
147 case '\v':
148 case '\f':
149 case '\r':
150 case ' ':
151 continue;
152 default:
153 goto break2;
154 }
155 break2:
156 if (*s == '0') {
157 #ifndef NO_HEX_FP /*{*/
158 {
159 static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
160 Long expt;
161 ULong bits[2];
162 switch(s[1]) {
163 case 'x':
164 case 'X':
165 {
166 #ifdef Honor_FLT_ROUNDS
167 FPI fpi1 = fpi;
168 fpi1.rounding = Rounding;
169 #else
170 #define fpi1 fpi
171 #endif
172 switch((i = gethex(&s, &fpi1, &expt, &bb, sign, loc)) & STRTOG_Retmask) {
173 case STRTOG_NoNumber:
174 s = s00;
175 sign = 0;
176 /* FALLTHROUGH */
177 case STRTOG_Zero:
178 break;
179 default:
180 if (bb) {
181 copybits(bits, fpi.nbits, bb);
182 Bfree(bb);
183 }
184 ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i);
185 }}
186 goto ret;
187 }
188 }
189 #endif /*}*/
190 nz0 = 1;
191 while(*++s == '0') ;
192 if (!*s)
193 goto ret;
194 }
195 s0 = s;
196 y = z = 0;
197 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
198 if (nd < 9)
199 y = 10*y + c - '0';
200 else if (nd < DBL_DIG + 2)
201 z = 10*z + c - '0';
202 nd0 = nd;
203 #ifdef USE_LOCALE
204 if (c == *decimalpoint) {
205 for(i = 1; decimalpoint[i]; ++i)
206 if (s[i] != decimalpoint[i])
207 goto dig_done;
208 s += i;
209 c = *s;
210 #else
211 if (c == '.') {
212 c = *++s;
213 #endif
214 #ifdef Honor_FLT_ROUNDS
215 decpt = 1;
216 #endif
217 if (!nd) {
218 for(; c == '0'; c = *++s)
219 nz++;
220 if (c > '0' && c <= '9') {
221 s0 = s;
222 nf += nz;
223 nz = 0;
224 goto have_dig;
225 }
226 goto dig_done;
227 }
228 for(; c >= '0' && c <= '9'; c = *++s) {
229 have_dig:
230 nz++;
231 if (c -= '0') {
232 nf += nz;
233 for(i = 1; i < nz; i++)
234 if (nd++ < 9)
235 y *= 10;
236 else if (nd <= DBL_DIG + 2)
237 z *= 10;
238 if (nd++ < 9)
239 y = 10*y + c;
240 else if (nd <= DBL_DIG + 2)
241 z = 10*z + c;
242 nz = 0;
243 }
244 }
245 }/*}*/
246 dig_done:
247 e = 0;
248 if (c == 'e' || c == 'E') {
249 if (!nd && !nz && !nz0) {
250 goto ret0;
251 }
252 s00 = s;
253 esign = 0;
254 switch(c = *++s) {
255 case '-':
256 esign = 1;
257 /* FALLTHROUGH */
258 case '+':
259 c = *++s;
260 }
261 if (c >= '0' && c <= '9') {
262 while(c == '0')
263 c = *++s;
264 if (c > '0' && c <= '9') {
265 L = c - '0';
266 s1 = s;
267 while((c = *++s) >= '0' && c <= '9')
268 L = 10*L + c - '0';
269 if (s - s1 > 8 || L > 19999)
270 /* Avoid confusion from exponents
271 * so large that e might overflow.
272 */
273 e = 19999; /* safe for 16 bit ints */
274 else
275 e = (int)L;
276 if (esign)
277 e = -e;
278 }
279 else
280 e = 0;
281 }
282 else
283 s = s00;
284 }
285 if (!nd) {
286 if (!nz && !nz0) {
287 #ifdef INFNAN_CHECK
288 /* Check for Nan and Infinity */
289 ULong bits[2];
290 static CONST FPI fpinan = /* only 52 explicit bits */
291 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
292 if (!decpt)
293 switch(c) {
294 case 'i':
295 case 'I':
296 if (match(&s,"nf")) {
297 --s;
298 if (!match(&s,"inity"))
299 ++s;
300 word0(&rv) = 0x7ff00000;
301 word1(&rv) = 0;
302 goto ret;
303 }
304 break;
305 case 'n':
306 case 'N':
307 if (match(&s, "an")) {
308 #ifndef No_Hex_NaN
309 if (*s == '(' /*)*/
310 && hexnan(&s, &fpinan, bits)
311 == STRTOG_NaNbits) {
312 word0(&rv) = 0x7ff00000 | bits[1];
313 word1(&rv) = bits[0];
314 }
315 else {
316 #endif
317 word0(&rv) = NAN_WORD0;
318 word1(&rv) = NAN_WORD1;
319 #ifndef No_Hex_NaN
320 }
321 #endif
322 goto ret;
323 }
324 }
325 #endif /* INFNAN_CHECK */
326 ret0:
327 s = s00;
328 sign = 0;
329 }
330 goto ret;
331 }
332 e1 = e -= nf;
333
334 /* Now we have nd0 digits, starting at s0, followed by a
335 * decimal point, followed by nd-nd0 digits. The number we're
336 * after is the integer represented by those digits times
337 * 10**e */
338
339 if (!nd0)
340 nd0 = nd;
341 k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
342 dval(&rv) = y;
343 if (k > 9) {
344 #ifdef SET_INEXACT
345 if (k > DBL_DIG)
346 oldinexact = get_inexact();
347 #endif
348 dval(&rv) = tens[k - 9] * dval(&rv) + z;
349 }
350 bd0 = 0;
351 if (nd <= DBL_DIG
352 #ifndef RND_PRODQUOT
353 #ifndef Honor_FLT_ROUNDS
354 /*CONSTCOND*/
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