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