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