]> git.saurik.com Git - apple/libc.git/blob - gdtoa/gdtoa-strtod-fbsd.c
projects / apple / libc.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Libc-583.tar.gz
[apple/libc.git] / gdtoa / gdtoa-strtod-fbsd.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 double aadj, aadj1, adj, rv, rv0;
77 Long L;
78 ULong y, z;
79 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
80 #ifdef SET_INEXACT
81 int inexact, oldinexact;
82 #endif
83 #ifdef USE_LOCALE
84 NORMALIZE_LOCALE(loc);
85 #ifdef NO_LOCALE_CACHE
86 char *decimalpoint = localeconv_l(loc)->decimal_point;
87 #else
88 char *decimalpoint;
89 static char *decimalpoint_cache;
90 if (!(s0 = decimalpoint_cache)) {
91 s0 = localeconv_l(loc)->decimal_point;
92 if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
93 strcpy(decimalpoint_cache, s0);
94 s0 = decimalpoint_cache;
95 }
96 }
97 decimalpoint = (char*)s0;
98 #endif
99 #endif
100 #ifdef Honor_FLT_ROUNDS /*{*/
101 int Rounding;
102 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
103 Rounding = Flt_Rounds;
104 #else /*}{*/
105 Rounding = 1;
106 switch(fegetround()) {
107 case FE_TOWARDZERO: Rounding = 0; break;
108 case FE_UPWARD: Rounding = 2; break;
109 case FE_DOWNWARD: Rounding = 3;
110 }
111 #endif /*}}*/
112 #endif /*}*/
113
114 sign = nz0 = nz = decpt = 0;
115 dval(rv) = 0.;
116 for(s = s00;;s++) switch(*s) {
117 case '-':
118 sign = 1;
119 /* no break */
120 case '+':
121 if (*++s)
122 goto break2;
123 /* no break */
124 case 0:
125 goto ret0;
126 case '\t':
127 case '\n':
128 case '\v':
129 case '\f':
130 case '\r':
131 case ' ':
132 continue;
133 default:
134 goto break2;
135 }
136 break2:
137 if (*s == '0') {
138 #ifndef NO_HEX_FP /*{*/
139 {
140 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
141 Long exp;
142 ULong bits[2];
143 switch(s[1]) {
144 case 'x':
145 case 'X':
146 {
147 #if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) /*{{*/
148 FPI fpi1 = fpi;
149 #ifdef Honor_FLT_ROUNDS /*{{*/
150 fpi1.rounding = Rounding;
151 #else /*}{*/
152 switch(fegetround()) {
153 case FE_TOWARDZERO: fpi1.rounding = 0; break;
154 case FE_UPWARD: fpi1.rounding = 2; break;
155 case FE_DOWNWARD: fpi1.rounding = 3;
156 }
157 #endif /*}}*/
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 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 e1 = e -= nf;
318
319 /* Now we have nd0 digits, starting at s0, followed by a
320 * decimal point, followed by nd-nd0 digits. The number we're
321 * after is the integer represented by those digits times
322 * 10**e */
323
324 if (!nd0)
325 nd0 = nd;
326 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
327 dval(rv) = y;
328 if (k > 9) {
329 #ifdef SET_INEXACT
330 if (k > DBL_DIG)
331 oldinexact = get_inexact();
332 #endif
333 dval(rv) = tens[k - 9] * dval(rv) + z;
334 }
335 bd0 = 0;
336 if (nd <= DBL_DIG
337 #ifndef RND_PRODQUOT
338 #ifndef Honor_FLT_ROUNDS
339 && Flt_Rounds == 1
340 #endif
341 #endif
342 ) {
343 if (!e)
344 goto ret;
345 if (e > 0) {
346 if (e <= Ten_pmax) {
347 #ifdef VAX
348 goto vax_ovfl_check;
349 #else
350 #ifdef Honor_FLT_ROUNDS
351 /* round correctly FLT_ROUNDS = 2 or 3 */
352 if (sign) {
353 rv = -rv;
354 sign = 0;
355 }
356 #endif
357 /* rv = */ rounded_product(dval(rv), tens[e]);
358 goto ret;
359 #endif
360 }
361 i = DBL_DIG - nd;
362 if (e <= Ten_pmax + i) {
363 /* A fancier test would sometimes let us do
364 * this for larger i values.
365 */
366 #ifdef Honor_FLT_ROUNDS
367 /* round correctly FLT_ROUNDS = 2 or 3 */
368 if (sign) {
369 rv = -rv;
370 sign = 0;
371 }
372 #endif
373 e -= i;
374 dval(rv) *= tens[i];
375 #ifdef VAX
376 /* VAX exponent range is so narrow we must
377 * worry about overflow here...
378 */
379 vax_ovfl_check:
380 word0(rv) -= P*Exp_msk1;
381 /* rv = */ rounded_product(dval(rv), tens[e]);
382 if ((word0(rv) & Exp_mask)
383 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
384 goto ovfl;
385 word0(rv) += P*Exp_msk1;
386 #else
387 /* rv = */ rounded_product(dval(rv), tens[e]);
388 #endif
389 goto ret;
390 }
391 }
392 #ifndef Inaccurate_Divide
393 else if (e >= -Ten_pmax) {
394 #ifdef Honor_FLT_ROUNDS
395 /* round correctly FLT_ROUNDS = 2 or 3 */
396 if (sign) {
397 rv = -rv;
398 sign = 0;
399 }
400 #endif
401 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
402 goto ret;
403 }
404 #endif
405 }
406 e1 += nd - k;
407
408 #ifdef IEEE_Arith
409 #ifdef SET_INEXACT
410 inexact = 1;
411 if (k <= DBL_DIG)
412 oldinexact = get_inexact();
413 #endif
414 #ifdef Avoid_Underflow
415 scale = 0;
416 #endif
417 #ifdef Honor_FLT_ROUNDS
418 if (Rounding >= 2) {
419 if (sign)
420 Rounding = Rounding == 2 ? 0 : 2;
421 else
422 if (Rounding != 2)
423 Rounding = 0;
424 }
425 #endif
426 #endif /*IEEE_Arith*/
427
428 /* Get starting approximation = rv * 10**e1 */
429
430 if (e1 > 0) {
431 if ( (i = e1 & 15) !=0)
432 dval(rv) *= tens[i];
433 if (e1 &= ~15) {
434 if (e1 > DBL_MAX_10_EXP) {
435 ovfl:
436 #ifndef NO_ERRNO
437 errno = ERANGE;
438 #endif
439 /* Can't trust HUGE_VAL */
440 #ifdef IEEE_Arith
441 #ifdef Honor_FLT_ROUNDS
442 switch(Rounding) {
443 case 0: /* toward 0 */
444 case 3: /* toward -infinity */
445 word0(rv) = Big0;
446 word1(rv) = Big1;
447 break;
448 default:
449 word0(rv) = Exp_mask;
450 word1(rv) = 0;
451 }
452 #else /*Honor_FLT_ROUNDS*/
453 word0(rv) = Exp_mask;
454 word1(rv) = 0;
455 #endif /*Honor_FLT_ROUNDS*/
456 #ifdef SET_INEXACT
457 /* set overflow bit */
458 dval(rv0) = 1e300;
459 dval(rv0) *= dval(rv0);
460 #endif
461 #else /*IEEE_Arith*/
462 word0(rv) = Big0;
463 word1(rv) = Big1;
464 #endif /*IEEE_Arith*/
465 if (bd0)
466 goto retfree;
467 goto ret;
468 }
469 e1 >>= 4;
470 for(j = 0; e1 > 1; j++, e1 >>= 1)
471 if (e1 & 1)
472 dval(rv) *= bigtens[j];
473 /* The last multiplication could overflow. */
474 word0(rv) -= P*Exp_msk1;
475 dval(rv) *= bigtens[j];
476 if ((z = word0(rv) & Exp_mask)
477 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
478 goto ovfl;
479 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
480 /* set to largest number */
481 /* (Can't trust DBL_MAX) */
482 word0(rv) = Big0;
483 word1(rv) = Big1;
484 }
485 else
486 word0(rv) += P*Exp_msk1;
487 }
488 }
489 else if (e1 < 0) {
490 e1 = -e1;
491 if ( (i = e1 & 15) !=0)
492 dval(rv) /= tens[i];
493 if (e1 >>= 4) {
494 if (e1 >= 1 << n_bigtens)
495 goto undfl;
496 #ifdef Avoid_Underflow
497 if (e1 & Scale_Bit)
498 scale = 2*P;
499 for(j = 0; e1 > 0; j++, e1 >>= 1)
500 if (e1 & 1)
501 dval(rv) *= tinytens[j];
502 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
503 >> Exp_shift)) > 0) {
504 /* scaled rv is denormal; zap j low bits */
505 if (j >= 32) {
506 word1(rv) = 0;
507 if (j >= 53)
508 word0(rv) = (P+2)*Exp_msk1;
509 else
510 word0(rv) &= 0xffffffff << j-32;
511 }
512 else
513 word1(rv) &= 0xffffffff << j;
514 }
515 #else
516 for(j = 0; e1 > 1; j++, e1 >>= 1)
517 if (e1 & 1)
518 dval(rv) *= tinytens[j];
519 /* The last multiplication could underflow. */
520 dval(rv0) = dval(rv);
521 dval(rv) *= tinytens[j];
522 if (!dval(rv)) {
523 dval(rv) = 2.*dval(rv0);
524 dval(rv) *= tinytens[j];
525 #endif
526 if (!dval(rv)) {
527 undfl:
528 dval(rv) = 0.;
529 #ifndef NO_ERRNO
530 errno = ERANGE;
531 #endif
532 if (bd0)
533 goto retfree;
534 goto ret;
535 }
536 #ifndef Avoid_Underflow
537 word0(rv) = Tiny0;
538 word1(rv) = Tiny1;
539 /* The refinement below will clean
540 * this approximation up.
541 */
542 }
543 #endif
544 }
545 }
546
547 /* Now the hard part -- adjusting rv to the correct value.*/
548
549 /* Put digits into bd: true value = bd * 10^e */
550
551 #ifdef USE_LOCALE
552 bd0 = s2b(s0, nd0, nd, y, strlen(decimalpoint));
553 #else
554 bd0 = s2b(s0, nd0, nd, y, 1);
555 #endif
556
557 for(;;) {
558 bd = Balloc(bd0->k);
559 Bcopy(bd, bd0);
560 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
561 bs = i2b(1);
562
563 if (e >= 0) {
564 bb2 = bb5 = 0;
565 bd2 = bd5 = e;
566 }
567 else {
568 bb2 = bb5 = -e;
569 bd2 = bd5 = 0;
570 }
571 if (bbe >= 0)
572 bb2 += bbe;
573 else
574 bd2 -= bbe;
575 bs2 = bb2;
576 #ifdef Honor_FLT_ROUNDS
577 if (Rounding != 1)
578 bs2++;
579 #endif
580 #ifdef Avoid_Underflow
581 j = bbe - scale;
582 i = j + bbbits - 1; /* logb(rv) */
583 if (i < Emin) /* denormal */
584 j += P - Emin;
585 else
586 j = P + 1 - bbbits;
587 #else /*Avoid_Underflow*/
588 #ifdef Sudden_Underflow
589 #ifdef IBM
590 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
591 #else
592 j = P + 1 - bbbits;
593 #endif
594 #else /*Sudden_Underflow*/
595 j = bbe;
596 i = j + bbbits - 1; /* logb(rv) */
597 if (i < Emin) /* denormal */
598 j += P - Emin;
599 else
600 j = P + 1 - bbbits;
601 #endif /*Sudden_Underflow*/
602 #endif /*Avoid_Underflow*/
603 bb2 += j;
604 bd2 += j;
605 #ifdef Avoid_Underflow
606 bd2 += scale;
607 #endif
608 i = bb2 < bd2 ? bb2 : bd2;
609 if (i > bs2)
610 i = bs2;
611 if (i > 0) {
612 bb2 -= i;
613 bd2 -= i;
614 bs2 -= i;
615 }
616 if (bb5 > 0) {
617 bs = pow5mult(bs, bb5);
618 bb1 = mult(bs, bb);
619 Bfree(bb);
620 bb = bb1;
621 }
622 if (bb2 > 0)
623 bb = lshift(bb, bb2);
624 if (bd5 > 0)
625 bd = pow5mult(bd, bd5);
626 if (bd2 > 0)
627 bd = lshift(bd, bd2);
628 if (bs2 > 0)
629 bs = lshift(bs, bs2);
630 delta = diff(bb, bd);
631 dsign = delta->sign;
632 delta->sign = 0;
633 i = cmp(delta, bs);
634 #ifdef Honor_FLT_ROUNDS
635 if (Rounding != 1) {
636 if (i < 0) {
637 /* Error is less than an ulp */
638 if (!delta->x[0] && delta->wds <= 1) {
639 /* exact */
640 #ifdef SET_INEXACT
641 inexact = 0;
642 #endif
643 break;
644 }
645 if (Rounding) {
646 if (dsign) {
647 adj = 1.;
648 goto apply_adj;
649 }
650 }
651 else if (!dsign) {
652 adj = -1.;
653 if (!word1(rv)
654 && !(word0(rv) & Frac_mask)) {
655 y = word0(rv) & Exp_mask;
656 #ifdef Avoid_Underflow
657 if (!scale || y > 2*P*Exp_msk1)
658 #else
659 if (y)
660 #endif
661 {
662 delta = lshift(delta,Log2P);
663 if (cmp(delta, bs) <= 0)
664 adj = -0.5;
665 }
666 }
667 apply_adj:
668 #ifdef Avoid_Underflow
669 if (scale && (y = word0(rv) & Exp_mask)
670 <= 2*P*Exp_msk1)
671 word0(adj) += (2*P+1)*Exp_msk1 - y;
672 #else
673 #ifdef Sudden_Underflow
674 if ((word0(rv) & Exp_mask) <=
675 P*Exp_msk1) {
676 word0(rv) += P*Exp_msk1;
677 dval(rv) += adj*ulp(dval(rv));
678 word0(rv) -= P*Exp_msk1;
679 }
680 else
681 #endif /*Sudden_Underflow*/
682 #endif /*Avoid_Underflow*/
683 dval(rv) += adj*ulp(dval(rv));
684 }
685 break;
686 }
687 adj = ratio(delta, bs);
688 if (adj < 1.)
689 adj = 1.;
690 if (adj <= 0x7ffffffe) {
691 /* adj = Rounding ? ceil(adj) : floor(adj); */
692 y = adj;
693 if (y != adj) {
694 if (!((Rounding>>1) ^ dsign))
695 y++;
696 adj = y;
697 }
698 }
699 #ifdef Avoid_Underflow
700 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
701 word0(adj) += (2*P+1)*Exp_msk1 - y;
702 #else
703 #ifdef Sudden_Underflow
704 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
705 word0(rv) += P*Exp_msk1;
706 adj *= ulp(dval(rv));
707 if (dsign)
708 dval(rv) += adj;
709 else
710 dval(rv) -= adj;
711 word0(rv) -= P*Exp_msk1;
712 goto cont;
713 }
714 #endif /*Sudden_Underflow*/
715 #endif /*Avoid_Underflow*/
716 adj *= ulp(dval(rv));
717 if (dsign) {
718 if (word0(rv) == Big0 && word1(rv) == Big1)
719 goto ovfl;
720 dval(rv) += adj;
721 }
722 else
723 dval(rv) -= adj;
724 goto cont;
725 }
726 #endif /*Honor_FLT_ROUNDS*/
727
728 if (i < 0) {
729 /* Error is less than half an ulp -- check for
730 * special case of mantissa a power of two.
731 */
732 if (dsign || word1(rv) || word0(rv) & Bndry_mask
733 #ifdef IEEE_Arith
734 #ifdef Avoid_Underflow
735 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
736 #else
737 || (word0(rv) & Exp_mask) <= Exp_msk1
738 #endif
739 #endif
740 ) {
741 #ifdef SET_INEXACT
742 if (!delta->x[0] && delta->wds <= 1)
743 inexact = 0;
744 #endif
745 break;
746 }
747 if (!delta->x[0] && delta->wds <= 1) {
748 /* exact result */
749 #ifdef SET_INEXACT
750 inexact = 0;
751 #endif
752 break;
753 }
754 delta = lshift(delta,Log2P);
755 if (cmp(delta, bs) > 0)
756 goto drop_down;
757 break;
758 }
759 if (i == 0) {
760 /* exactly half-way between */
761 if (dsign) {
762 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
763 && word1(rv) == (
764 #ifdef Avoid_Underflow
765 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
766 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
767 #endif
768 0xffffffff)) {
769 /*boundary case -- increment exponent*/
770 word0(rv) = (word0(rv) & Exp_mask)
771 + Exp_msk1
772 #ifdef IBM
773 | Exp_msk1 >> 4
774 #endif
775 ;
776 word1(rv) = 0;
777 #ifdef Avoid_Underflow
778 dsign = 0;
779 #endif
780 break;
781 }
782 }
783 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
784 drop_down:
785 /* boundary case -- decrement exponent */
786 #ifdef Sudden_Underflow /*{{*/
787 L = word0(rv) & Exp_mask;
788 #ifdef IBM
789 if (L < Exp_msk1)
790 #else
791 #ifdef Avoid_Underflow
792 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
793 #else
794 if (L <= Exp_msk1)
795 #endif /*Avoid_Underflow*/
796 #endif /*IBM*/
797 goto undfl;
798 L -= Exp_msk1;
799 #else /*Sudden_Underflow}{*/
800 #ifdef Avoid_Underflow
801 if (scale) {
802 L = word0(rv) & Exp_mask;
803 if (L <= (2*P+1)*Exp_msk1) {
804 if (L > (P+2)*Exp_msk1)
805 /* round even ==> */
806 /* accept rv */
807 break;
808 /* rv = smallest denormal */
809 goto undfl;
810 }
811 }
812 #endif /*Avoid_Underflow*/
813 L = (word0(rv) & Exp_mask) - Exp_msk1;
814 #endif /*Sudden_Underflow}}*/
815 word0(rv) = L | Bndry_mask1;
816 word1(rv) = 0xffffffff;
817 #ifdef IBM
818 goto cont;
819 #else
820 break;
821 #endif
822 }
823 #ifndef ROUND_BIASED
824 if (!(word1(rv) & LSB))
825 break;
826 #endif
827 if (dsign)
828 dval(rv) += ulp(dval(rv));
829 #ifndef ROUND_BIASED
830 else {
831 dval(rv) -= ulp(dval(rv));
832 #ifndef Sudden_Underflow
833 if (!dval(rv))
834 goto undfl;
835 #endif
836 }
837 #ifdef Avoid_Underflow
838 dsign = 1 - dsign;
839 #endif
840 #endif
841 break;
842 }
843 if ((aadj = ratio(delta, bs)) <= 2.) {
844 if (dsign)
845 aadj = aadj1 = 1.;
846 else if (word1(rv) || word0(rv) & Bndry_mask) {
847 #ifndef Sudden_Underflow
848 if (word1(rv) == Tiny1 && !word0(rv))
849 goto undfl;
850 #endif
851 aadj = 1.;
852 aadj1 = -1.;
853 }
854 else {
855 /* special case -- power of FLT_RADIX to be */
856 /* rounded down... */
857
858 if (aadj < 2./FLT_RADIX)
859 aadj = 1./FLT_RADIX;
860 else
861 aadj *= 0.5;
862 aadj1 = -aadj;
863 }
864 }
865 else {
866 aadj *= 0.5;
867 aadj1 = dsign ? aadj : -aadj;
868 #ifdef Check_FLT_ROUNDS
869 switch(Rounding) {
870 case 2: /* towards +infinity */
871 aadj1 -= 0.5;
872 break;
873 case 0: /* towards 0 */
874 case 3: /* towards -infinity */
875 aadj1 += 0.5;
876 }
877 #else
878 if (Flt_Rounds == 0)
879 aadj1 += 0.5;
880 #endif /*Check_FLT_ROUNDS*/
881 }
882 y = word0(rv) & Exp_mask;
883
884 /* Check for overflow */
885
886 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
887 dval(rv0) = dval(rv);
888 word0(rv) -= P*Exp_msk1;
889 adj = aadj1 * ulp(dval(rv));
890 dval(rv) += adj;
891 if ((word0(rv) & Exp_mask) >=
892 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
893 if (word0(rv0) == Big0 && word1(rv0) == Big1)
894 goto ovfl;
895 word0(rv) = Big0;
896 word1(rv) = Big1;
897 goto cont;
898 }
899 else
900 word0(rv) += P*Exp_msk1;
901 }
902 else {
903 #ifdef Avoid_Underflow
904 if (scale && y <= 2*P*Exp_msk1) {
905 if (aadj <= 0x7fffffff) {
906 if ((z = aadj) <= 0)
907 z = 1;
908 aadj = z;
909 aadj1 = dsign ? aadj : -aadj;
910 }
911 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
912 }
913 adj = aadj1 * ulp(dval(rv));
914 dval(rv) += adj;
915 #else
916 #ifdef Sudden_Underflow
917 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
918 dval(rv0) = dval(rv);
919 word0(rv) += P*Exp_msk1;
920 adj = aadj1 * ulp(dval(rv));
921 dval(rv) += adj;
922 #ifdef IBM
923 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
924 #else
925 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
926 #endif
927 {
928 if (word0(rv0) == Tiny0
929 && word1(rv0) == Tiny1)
930 goto undfl;
931 word0(rv) = Tiny0;
932 word1(rv) = Tiny1;
933 goto cont;
934 }
935 else
936 word0(rv) -= P*Exp_msk1;
937 }
938 else {
939 adj = aadj1 * ulp(dval(rv));
940 dval(rv) += adj;
941 }
942 #else /*Sudden_Underflow*/
943 /* Compute adj so that the IEEE rounding rules will
944 * correctly round rv + adj in some half-way cases.
945 * If rv * ulp(rv) is denormalized (i.e.,
946 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
947 * trouble from bits lost to denormalization;
948 * example: 1.2e-307 .
949 */
950 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
951 aadj1 = (double)(int)(aadj + 0.5);
952 if (!dsign)
953 aadj1 = -aadj1;
954 }
955 adj = aadj1 * ulp(dval(rv));
956 dval(rv) += adj;
957 #endif /*Sudden_Underflow*/
958 #endif /*Avoid_Underflow*/
959 }
960 z = word0(rv) & Exp_mask;
961 #ifndef SET_INEXACT
962 #ifdef Avoid_Underflow
963 if (!scale)
964 #endif
965 if (y == z) {
966 /* Can we stop now? */
967 L = (Long)aadj;
968 aadj -= L;
969 /* The tolerances below are conservative. */
970 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
971 if (aadj < .4999999 || aadj > .5000001)
972 break;
973 }
974 else if (aadj < .4999999/FLT_RADIX)
975 break;
976 }
977 #endif
978 cont:
979 Bfree(bb);
980 Bfree(bd);
981 Bfree(bs);
982 Bfree(delta);
983 }
984 #ifdef SET_INEXACT
985 if (inexact) {
986 if (!oldinexact) {
987 word0(rv0) = Exp_1 + (70 << Exp_shift);
988 word1(rv0) = 0;
989 dval(rv0) += 1.;
990 }
991 }
992 else if (!oldinexact)
993 clear_inexact();
994 #endif
995 #ifdef Avoid_Underflow
996 if (scale) {
997 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
998 word1(rv0) = 0;
999 dval(rv) *= dval(rv0);
1000 #ifndef NO_ERRNO
1001 /* try to avoid the bug of testing an 8087 register value */
1002 #if defined(IEEE_Arith) && __DARWIN_UNIX03
1003 if (!(word0(rv) & Exp_mask))
1004 #else
1005 if (word0(rv) == 0 && word1(rv) == 0)
1006 #endif
1007 errno = ERANGE;
1008 #endif
1009 }
1010 #endif /* Avoid_Underflow */
1011 #ifdef SET_INEXACT
1012 if (inexact && !(word0(rv) & Exp_mask)) {
1013 /* set underflow bit */
1014 dval(rv0) = 1e-300;
1015 dval(rv0) *= dval(rv0);
1016 }
1017 #endif
1018 retfree:
1019 Bfree(bb);
1020 Bfree(bd);
1021 Bfree(bs);
1022 Bfree(bd0);
1023 Bfree(delta);
1024 ret:
1025 if (se)
1026 *se = (char *)s;
1027 return sign ? -dval(rv) : dval(rv);
1028 }
1029
1030 double
1031 strtod
1032 #ifdef KR_headers
1033 (s00, se) CONST char *s00; char **se;
1034 #else
1035 (CONST char *s00, char **se)
1036 #endif
1037 {
1038 return strtod_l(s00, se, __current_locale());
1039 }
mirror of libc