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