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 "gdtoaimp.h"
33
34 #ifdef USE_LOCALE
35 #include "locale.h"
36 #endif
37
38 static CONST int
39 fivesbits[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21,
40 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
41 47, 49, 52
42 #ifdef VAX
43 , 54, 56
44 #endif
45 };
46
47 Bigint *
48 #ifdef KR_headers
increment(b)49 increment(b) Bigint *b;
50 #else
51 increment(Bigint *b)
52 #endif
53 {
54 ULong *x, *xe;
55 Bigint *b1;
56 #ifdef Pack_16
57 ULong carry = 1, y;
58 #endif
59
60 x = b->x;
61 xe = x + b->wds;
62 #ifdef Pack_32
63 do {
64 if (*x < (ULong)0xffffffffL) {
65 ++*x;
66 return b;
67 }
68 *x++ = 0;
69 } while(x < xe);
70 #else
71 do {
72 y = *x + carry;
73 carry = y >> 16;
74 *x++ = y & 0xffff;
75 if (!carry)
76 return b;
77 } while(x < xe);
78 if (carry)
79 #endif
80 {
81 if (b->wds >= b->maxwds) {
82 b1 = Balloc(b->k+1);
83 Bcopy(b1,b);
84 Bfree(b);
85 b = b1;
86 }
87 b->x[b->wds++] = 1;
88 }
89 return b;
90 }
91
92 void
93 #ifdef KR_headers
decrement(b)94 decrement(b) Bigint *b;
95 #else
96 decrement(Bigint *b)
97 #endif
98 {
99 ULong *x, *xe;
100 #ifdef Pack_16
101 ULong borrow = 1, y;
102 #endif
103
104 x = b->x;
105 xe = x + b->wds;
106 #ifdef Pack_32
107 do {
108 if (*x) {
109 --*x;
110 break;
111 }
112 *x++ = 0xffffffffL;
113 }
114 while(x < xe);
115 #else
116 do {
117 y = *x - borrow;
118 borrow = (y & 0x10000) >> 16;
119 *x++ = y & 0xffff;
120 } while(borrow && x < xe);
121 #endif
122 }
123
124 static int
125 #ifdef KR_headers
all_on(b,n)126 all_on(b, n) Bigint *b; int n;
127 #else
128 all_on(Bigint *b, int n)
129 #endif
130 {
131 ULong *x, *xe;
132
133 x = b->x;
134 xe = x + (n >> kshift);
135 while(x < xe)
136 if ((*x++ & ALL_ON) != ALL_ON)
137 return 0;
138 if (n &= kmask)
139 return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON;
140 return 1;
141 }
142
143 Bigint *
144 #ifdef KR_headers
set_ones(b,n)145 set_ones(b, n) Bigint *b; int n;
146 #else
147 set_ones(Bigint *b, int n)
148 #endif
149 {
150 int k;
151 ULong *x, *xe;
152
153 k = (n + ((1 << kshift) - 1)) >> kshift;
154 if (b->k < k) {
155 Bfree(b);
156 b = Balloc(k);
157 }
158 k = n >> kshift;
159 if (n &= kmask)
160 k++;
161 b->wds = k;
162 x = b->x;
163 xe = x + k;
164 while(x < xe)
165 *x++ = ALL_ON;
166 if (n)
167 x[-1] >>= ULbits - n;
168 return b;
169 }
170
171 static int
rvOK(d,fpi,exp,bits,exact,rd,irv)172 rvOK
173 #ifdef KR_headers
174 (d, fpi, exp, bits, exact, rd, irv)
175 U *d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
176 #else
177 (U *d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
178 #endif
179 {
180 Bigint *b;
181 ULong carry, inex, lostbits;
182 int bdif, e, j, k, k1, nb, rv;
183
184 carry = rv = 0;
185 b = d2b(dval(d), &e, &bdif);
186 bdif -= nb = fpi->nbits;
187 e += bdif;
188 if (bdif <= 0) {
189 if (exact)
190 goto trunc;
191 goto ret;
192 }
193 if (P == nb) {
194 if (
195 #ifndef IMPRECISE_INEXACT
196 exact &&
197 #endif
198 fpi->rounding ==
199 #ifdef RND_PRODQUOT
200 FPI_Round_near
201 #else
202 Flt_Rounds
203 #endif
204 ) goto trunc;
205 goto ret;
206 }
207 switch(rd) {
208 case 1: /* round down (toward -Infinity) */
209 goto trunc;
210 case 2: /* round up (toward +Infinity) */
211 break;
212 default: /* round near */
213 k = bdif - 1;
214 if (k < 0)
215 goto trunc;
216 if (!k) {
217 if (!exact)
218 goto ret;
219 if (b->x[0] & 2)
220 break;
221 goto trunc;
222 }
223 if (b->x[k>>kshift] & ((ULong)1 << (k & kmask)))
224 break;
225 goto trunc;
226 }
227 /* "break" cases: round up 1 bit, then truncate; bdif > 0 */
228 carry = 1;
229 trunc:
230 inex = lostbits = 0;
231 if (bdif > 0) {
232 if ( (lostbits = any_on(b, bdif)) !=0)
233 inex = STRTOG_Inexlo;
234 rshift(b, bdif);
235 if (carry) {
236 inex = STRTOG_Inexhi;
237 b = increment(b);
238 if ( (j = nb & kmask) !=0)
239 j = ULbits - j;
240 if (hi0bits(b->x[b->wds - 1]) != j) {
241 if (!lostbits)
242 lostbits = b->x[0] & 1;
243 rshift(b, 1);
244 e++;
245 }
246 }
247 }
248 else if (bdif < 0)
249 b = lshift(b, -bdif);
250 if (e < fpi->emin) {
251 k = fpi->emin - e;
252 e = fpi->emin;
253 if (k > nb || fpi->sudden_underflow) {
254 b->wds = inex = 0;
255 *irv = STRTOG_Underflow | STRTOG_Inexlo;
256 }
257 else {
258 k1 = k - 1;
259 if (k1 > 0 && !lostbits)
260 lostbits = any_on(b, k1);
261 if (!lostbits && !exact)
262 goto ret;
263 lostbits |=
264 carry = b->x[k1>>kshift] & (1 << (k1 & kmask));
265 rshift(b, k);
266 *irv = STRTOG_Denormal;
267 if (carry) {
268 b = increment(b);
269 inex = STRTOG_Inexhi | STRTOG_Underflow;
270 }
271 else if (lostbits)
272 inex = STRTOG_Inexlo | STRTOG_Underflow;
273 }
274 }
275 else if (e > fpi->emax) {
276 e = fpi->emax + 1;
277 *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
278 #ifndef NO_ERRNO
279 errno = ERANGE;
280 #endif
281 b->wds = inex = 0;
282 }
283 *exp = e;
284 copybits(bits, nb, b);
285 *irv |= inex;
286 rv = 1;
287 ret:
288 Bfree(b);
289 return rv;
290 }
291
292 static int
293 #ifdef KR_headers
mantbits(d)294 mantbits(d) U *d;
295 #else
296 mantbits(U *d)
297 #endif
298 {
299 ULong L;
300 #ifdef VAX
301 L = word1(d) << 16 | word1(d) >> 16;
302 if (L)
303 #else
304 if ( (L = word1(d)) !=0)
305 #endif
306 return P - lo0bits(&L);
307 #ifdef VAX
308 L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;
309 #else
310 L = word0(d) | Exp_msk1;
311 #endif
312 return P - 32 - lo0bits(&L);
313 }
314
315 int
strtodg_l(s00,se,fpi,exp,bits,loc)316 strtodg_l
317 #ifdef KR_headers
318 (s00, se, fpi, exp, bits, loc)
319 CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits; locale_t loc;
320 #else
321 (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits, locale_t loc)
322 #endif
323 {
324 int abe, abits, asub;
325 int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
326 int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
327 int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
328 int sudden_underflow;
329 CONST char *s, *s0, *s1;
330 double adj0, tol;
331 Long L;
332 U adj, rv;
333 ULong *b, *be, y, z;
334 Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
335 #ifdef USE_LOCALE /*{{*/
336 #ifdef NO_LOCALE_CACHE
337 char *decimalpoint = localeconv_l(loc)->decimal_point;
338 int dplen = strlen(decimalpoint);
339 #else
340 char *decimalpoint;
341 static char *decimalpoint_cache;
342 static int dplen;
343 if (!(s0 = decimalpoint_cache)) {
344 s0 = localeconv_l(loc)->decimal_point;
345 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
346 strcpy(decimalpoint_cache, s0);
347 s0 = decimalpoint_cache;
348 }
349 dplen = strlen(s0);
350 }
351 decimalpoint = (char*)s0;
352 #endif /*NO_LOCALE_CACHE*/
353 #else /*USE_LOCALE}{*/
354 #define dplen 1
355 #endif /*USE_LOCALE}}*/
356
357 irv = STRTOG_Zero;
358 denorm = sign = nz0 = nz = 0;
359 dval(&rv) = 0.;
360 rvb = 0;
361 nbits = fpi->nbits;
362 for(s = s00;;s++) switch(*s) {
363 case '-':
364 sign = 1;
365 /* no break */
366 case '+':
367 if (*++s)
368 goto break2;
369 /* no break */
370 case 0:
371 sign = 0;
372 irv = STRTOG_NoNumber;
373 s = s00;
374 goto ret;
375 case '\t':
376 case '\n':
377 case '\v':
378 case '\f':
379 case '\r':
380 case ' ':
381 continue;
382 default:
383 goto break2;
384 }
385 break2:
386 if (*s == '0') {
387 #ifndef NO_HEX_FP
388 switch(s[1]) {
389 case 'x':
390 case 'X':
391 irv = gethex(&s, fpi, exp, &rvb, sign);
392 if (irv == STRTOG_NoNumber) {
393 s = s00;
394 sign = 0;
395 }
396 goto ret;
397 }
398 #endif
399 nz0 = 1;
400 while(*++s == '0') ;
401 if (!*s)
402 goto ret;
403 }
404 sudden_underflow = fpi->sudden_underflow;
405 s0 = s;
406 y = z = 0;
407 for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
408 if (nd < 9)
409 y = 10*y + c - '0';
410 else if (nd < 16)
411 z = 10*z + c - '0';
412 nd0 = nd;
413 #ifdef USE_LOCALE
414 if (c == *decimalpoint) {
415 for(i = 1; decimalpoint[i]; ++i)
416 if (s[i] != decimalpoint[i])
417 goto dig_done;
418 s += i;
419 c = *s;
420 #else
421 if (c == '.') {
422 c = *++s;
423 #endif
424 decpt = 1;
425 if (!nd) {
426 for(; c == '0'; c = *++s)
427 nz++;
428 if (c > '0' && c <= '9') {
429 s0 = s;
430 nf += nz;
431 nz = 0;
432 goto have_dig;
433 }
434 goto dig_done;
435 }
436 for(; c >= '0' && c <= '9'; c = *++s) {
437 have_dig:
438 nz++;
439 if (c -= '0') {
440 nf += nz;
441 for(i = 1; i < nz; i++)
442 if (nd++ < 9)
443 y *= 10;
444 else if (nd <= DBL_DIG + 1)
445 z *= 10;
446 if (nd++ < 9)
447 y = 10*y + c;
448 else if (nd <= DBL_DIG + 1)
449 z = 10*z + c;
450 nz = 0;
451 }
452 }
453 }/*}*/
454 dig_done:
455 e = 0;
456 if (c == 'e' || c == 'E') {
457 if (!nd && !nz && !nz0) {
458 irv = STRTOG_NoNumber;
459 s = s00;
460 goto ret;
461 }
462 s00 = s;
463 esign = 0;
464 switch(c = *++s) {
465 case '-':
466 esign = 1;
467 case '+':
468 c = *++s;
469 }
470 if (c >= '0' && c <= '9') {
471 while(c == '0')
472 c = *++s;
473 if (c > '0' && c <= '9') {
474 L = c - '0';
475 s1 = s;
476 while((c = *++s) >= '0' && c <= '9')
477 L = 10*L + c - '0';
478 if (s - s1 > 8 || L > 19999)
479 /* Avoid confusion from exponents
480 * so large that e might overflow.
481 */
482 e = 19999; /* safe for 16 bit ints */
483 else
484 e = (int)L;
485 if (esign)
486 e = -e;
487 }
488 else
489 e = 0;
490 }
491 else
492 s = s00;
493 }
494 if (!nd) {
495 if (!nz && !nz0) {
496 #ifdef INFNAN_CHECK
497 /* Check for Nan and Infinity */
498 if (!decpt)
499 switch(c) {
500 case 'i':
501 case 'I':
502 if (match(&s,"nf")) {
503 --s;
504 if (!match(&s,"inity"))
505 ++s;
506 irv = STRTOG_Infinite;
507 goto infnanexp;
508 }
509 break;
510 case 'n':
511 case 'N':
512 if (match(&s, "an")) {
513 irv = STRTOG_NaN;
514 *exp = fpi->emax + 1;
515 #ifndef No_Hex_NaN
516 if (*s == '(') /*)*/
517 irv = hexnan(&s, fpi, bits);
518 #endif
519 goto infnanexp;
520 }
521 }
522 #endif /* INFNAN_CHECK */
523 irv = STRTOG_NoNumber;
524 s = s00;
525 }
526 goto ret;
527 }
528
529 irv = STRTOG_Normal;
530 e1 = e -= nf;
531 rd = 0;
532 switch(fpi->rounding & 3) {
533 case FPI_Round_up:
534 rd = 2 - sign;
535 break;
536 case FPI_Round_zero:
537 rd = 1;
538 break;
539 case FPI_Round_down:
540 rd = 1 + sign;
541 }
542
543 /* Now we have nd0 digits, starting at s0, followed by a
544 * decimal point, followed by nd-nd0 digits. The number we're
545 * after is the integer represented by those digits times
546 * 10**e */
547
548 if (!nd0)
549 nd0 = nd;
550 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
551 dval(&rv) = y;
552 if (k > 9)
553 dval(&rv) = tens[k - 9] * dval(&rv) + z;
554 bd0 = 0;
555 if (nbits <= P && nd <= DBL_DIG) {
556 if (!e) {
557 if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv))
558 goto ret;
559 }
560 else if (e > 0) {
561 if (e <= Ten_pmax) {
562 #ifdef VAX
563 goto vax_ovfl_check;
564 #else
565 i = fivesbits[e] + mantbits(&rv) <= P;
566 /* rv = */ rounded_product(dval(&rv), tens[e]);
567 if (rvOK(&rv, fpi, exp, bits, i, rd, &irv))
568 goto ret;
569 e1 -= e;
570 goto rv_notOK;
571 #endif
572 }
573 i = DBL_DIG - nd;
574 if (e <= Ten_pmax + i) {
575 /* A fancier test would sometimes let us do
576 * this for larger i values.
577 */
578 e2 = e - i;
579 e1 -= i;
580 dval(&rv) *= tens[i];
581 #ifdef VAX
582 /* VAX exponent range is so narrow we must
583 * worry about overflow here...
584 */
585 vax_ovfl_check:
586 dval(&adj) = dval(&rv);
587 word0(&adj) -= P*Exp_msk1;
588 /* adj = */ rounded_product(dval(&adj), tens[e2]);
589 if ((word0(&adj) & Exp_mask)
590 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
591 goto rv_notOK;
592 word0(&adj) += P*Exp_msk1;
593 dval(&rv) = dval(&adj);
594 #else
595 /* rv = */ rounded_product(dval(&rv), tens[e2]);
596 #endif
597 if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
598 goto ret;
599 e1 -= e2;
600 }
601 }
602 #ifndef Inaccurate_Divide
603 else if (e >= -Ten_pmax) {
604 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
605 if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
606 goto ret;
607 e1 -= e;
608 }
609 #endif
610 }
611 rv_notOK:
612 e1 += nd - k;
613
614 /* Get starting approximation = rv * 10**e1 */
615
616 e2 = 0;
617 if (e1 > 0) {
618 if ( (i = e1 & 15) !=0)
619 dval(&rv) *= tens[i];
620 if (e1 &= ~15) {
621 e1 >>= 4;
622 while(e1 >= (1 << (n_bigtens-1))) {
623 e2 += ((word0(&rv) & Exp_mask)
624 >> Exp_shift1) - Bias;
625 word0(&rv) &= ~Exp_mask;
626 word0(&rv) |= Bias << Exp_shift1;
627 dval(&rv) *= bigtens[n_bigtens-1];
628 e1 -= 1 << (n_bigtens-1);
629 }
630 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
631 word0(&rv) &= ~Exp_mask;
632 word0(&rv) |= Bias << Exp_shift1;
633 for(j = 0; e1 > 0; j++, e1 >>= 1)
634 if (e1 & 1)
635 dval(&rv) *= bigtens[j];
636 }
637 }
638 else if (e1 < 0) {
639 e1 = -e1;
640 if ( (i = e1 & 15) !=0)
641 dval(&rv) /= tens[i];
642 if (e1 &= ~15) {
643 e1 >>= 4;
644 while(e1 >= (1 << (n_bigtens-1))) {
645 e2 += ((word0(&rv) & Exp_mask)
646 >> Exp_shift1) - Bias;
647 word0(&rv) &= ~Exp_mask;
648 word0(&rv) |= Bias << Exp_shift1;
649 dval(&rv) *= tinytens[n_bigtens-1];
650 e1 -= 1 << (n_bigtens-1);
651 }
652 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
653 word0(&rv) &= ~Exp_mask;
654 word0(&rv) |= Bias << Exp_shift1;
655 for(j = 0; e1 > 0; j++, e1 >>= 1)
656 if (e1 & 1)
657 dval(&rv) *= tinytens[j];
658 }
659 }
660 #ifdef IBM
661 /* e2 is a correction to the (base 2) exponent of the return
662 * value, reflecting adjustments above to avoid overflow in the
663 * native arithmetic. For native IBM (base 16) arithmetic, we
664 * must multiply e2 by 4 to change from base 16 to 2.
665 */
666 e2 <<= 2;
667 #endif
668 rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */
669 rve += e2;
670 if ((j = rvbits - nbits) > 0) {
671 rshift(rvb, j);
672 rvbits = nbits;
673 rve += j;
674 }
675 bb0 = 0; /* trailing zero bits in rvb */
676 e2 = rve + rvbits - nbits;
677 if (e2 > fpi->emax + 1)
678 goto huge;
679 rve1 = rve + rvbits - nbits;
680 if (e2 < (emin = fpi->emin)) {
681 denorm = 1;
682 j = rve - emin;
683 if (j > 0) {
684 rvb = lshift(rvb, j);
685 rvbits += j;
686 }
687 else if (j < 0) {
688 rvbits += j;
689 if (rvbits <= 0) {
690 if (rvbits < -1) {
691 ufl:
692 rvb->wds = 0;
693 rvb->x[0] = 0;
694 *exp = emin;
695 irv = STRTOG_Underflow | STRTOG_Inexlo;
696 goto ret;
697 }
698 rvb->x[0] = rvb->wds = rvbits = 1;
699 }
700 else
701 rshift(rvb, -j);
702 }
703 rve = rve1 = emin;
704 if (sudden_underflow && e2 + 1 < emin)
705 goto ufl;
706 }
707
708 /* Now the hard part -- adjusting rv to the correct value.*/
709
710 /* Put digits into bd: true value = bd * 10^e */
711
712 bd0 = s2b(s0, nd0, nd, y, dplen);
713
714 for(;;) {
715 bd = Balloc(bd0->k);
716 Bcopy(bd, bd0);
717 bb = Balloc(rvb->k);
718 Bcopy(bb, rvb);
719 bbbits = rvbits - bb0;
720 bbe = rve + bb0;
721 bs = i2b(1);
722
723 if (e >= 0) {
724 bb2 = bb5 = 0;
725 bd2 = bd5 = e;
726 }
727 else {
728 bb2 = bb5 = -e;
729 bd2 = bd5 = 0;
730 }
731 if (bbe >= 0)
732 bb2 += bbe;
733 else
734 bd2 -= bbe;
735 bs2 = bb2;
736 j = nbits + 1 - bbbits;
737 i = bbe + bbbits - nbits;
738 if (i < emin) /* denormal */
739 j += i - emin;
740 bb2 += j;
741 bd2 += j;
742 i = bb2 < bd2 ? bb2 : bd2;
743 if (i > bs2)
744 i = bs2;
745 if (i > 0) {
746 bb2 -= i;
747 bd2 -= i;
748 bs2 -= i;
749 }
750 if (bb5 > 0) {
751 bs = pow5mult(bs, bb5);
752 bb1 = mult(bs, bb);
753 Bfree(bb);
754 bb = bb1;
755 }
756 bb2 -= bb0;
757 if (bb2 > 0)
758 bb = lshift(bb, bb2);
759 else if (bb2 < 0)
760 rshift(bb, -bb2);
761 if (bd5 > 0)
762 bd = pow5mult(bd, bd5);
763 if (bd2 > 0)
764 bd = lshift(bd, bd2);
765 if (bs2 > 0)
766 bs = lshift(bs, bs2);
767 asub = 1;
768 inex = STRTOG_Inexhi;
769 delta = diff(bb, bd);
770 if (delta->wds <= 1 && !delta->x[0])
771 break;
772 dsign = delta->sign;
773 delta->sign = finished = 0;
774 L = 0;
775 i = cmp(delta, bs);
776 if (rd && i <= 0) {
777 irv = STRTOG_Normal;
778 if ( (finished = dsign ^ (rd&1)) !=0) {
779 if (dsign != 0) {
780 irv |= STRTOG_Inexhi;
781 goto adj1;
782 }
783 irv |= STRTOG_Inexlo;
784 if (rve1 == emin)
785 goto adj1;
786 for(i = 0, j = nbits; j >= ULbits;
787 i++, j -= ULbits) {
788 if (rvb->x[i] & ALL_ON)
789 goto adj1;
790 }
791 if (j > 1 && lo0bits(rvb->x + i) < j - 1)
792 goto adj1;
793 rve = rve1 - 1;
794 rvb = set_ones(rvb, rvbits = nbits);
795 break;
796 }
797 irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
798 break;
799 }
800 if (i < 0) {
801 /* Error is less than half an ulp -- check for
802 * special case of mantissa a power of two.
803 */
804 irv = dsign
805 ? STRTOG_Normal | STRTOG_Inexlo
806 : STRTOG_Normal | STRTOG_Inexhi;
807 if (dsign || bbbits > 1 || denorm || rve1 == emin)
808 break;
809 delta = lshift(delta,1);
810 if (cmp(delta, bs) > 0) {
811 irv = STRTOG_Normal | STRTOG_Inexlo;
812 goto drop_down;
813 }
814 break;
815 }
816 if (i == 0) {
817 /* exactly half-way between */
818 if (dsign) {
819 if (denorm && all_on(rvb, rvbits)) {
820 /*boundary case -- increment exponent*/
821 rvb->wds = 1;
822 rvb->x[0] = 1;
823 rve = emin + nbits - (rvbits = 1);
824 irv = STRTOG_Normal | STRTOG_Inexhi;
825 denorm = 0;
826 break;
827 }
828 irv = STRTOG_Normal | STRTOG_Inexlo;
829 }
830 else if (bbbits == 1) {
831 irv = STRTOG_Normal;
832 drop_down:
833 /* boundary case -- decrement exponent */
834 if (rve1 == emin) {
835 irv = STRTOG_Normal | STRTOG_Inexhi;
836 if (rvb->wds == 1 && rvb->x[0] == 1)
837 sudden_underflow = 1;
838 break;
839 }
840 rve -= nbits;
841 rvb = set_ones(rvb, rvbits = nbits);
842 break;
843 }
844 else
845 irv = STRTOG_Normal | STRTOG_Inexhi;
846 if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
847 break;
848 if (dsign) {
849 rvb = increment(rvb);
850 j = kmask & (ULbits - (rvbits & kmask));
851 if (hi0bits(rvb->x[rvb->wds - 1]) != j)
852 rvbits++;
853 irv = STRTOG_Normal | STRTOG_Inexhi;
854 }
855 else {
856 if (bbbits == 1)
857 goto undfl;
858 decrement(rvb);
859 irv = STRTOG_Normal | STRTOG_Inexlo;
860 }
861 break;
862 }
863 if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
864 adj1:
865 inex = STRTOG_Inexlo;
866 if (dsign) {
867 asub = 0;
868 inex = STRTOG_Inexhi;
869 }
870 else if (denorm && bbbits <= 1) {
871 undfl:
872 rvb->wds = 0;
873 rve = emin;
874 irv = STRTOG_Underflow | STRTOG_Inexlo;
875 break;
876 }
877 adj0 = dval(&adj) = 1.;
878 }
879 else {
880 adj0 = dval(&adj) *= 0.5;
881 if (dsign) {
882 asub = 0;
883 inex = STRTOG_Inexlo;
884 }
885 if (dval(&adj) < 2147483647.) {
886 L = adj0;
887 adj0 -= L;
888 switch(rd) {
889 case 0:
890 if (adj0 >= .5)
891 goto inc_L;
892 break;
893 case 1:
894 if (asub && adj0 > 0.)
895 goto inc_L;
896 break;
897 case 2:
898 if (!asub && adj0 > 0.) {
899 inc_L:
900 L++;
901 inex = STRTOG_Inexact - inex;
902 }
903 }
904 dval(&adj) = L;
905 }
906 }
907 y = rve + rvbits;
908
909 /* adj *= ulp(dval(&rv)); */
910 /* if (asub) rv -= adj; else rv += adj; */
911
912 if (!denorm && rvbits < nbits) {
913 rvb = lshift(rvb, j = nbits - rvbits);
914 rve -= j;
915 rvbits = nbits;
916 }
917 ab = d2b(dval(&adj), &abe, &abits);
918 if (abe < 0)
919 rshift(ab, -abe);
920 else if (abe > 0)
921 ab = lshift(ab, abe);
922 rvb0 = rvb;
923 if (asub) {
924 /* rv -= adj; */
925 j = hi0bits(rvb->x[rvb->wds-1]);
926 rvb = diff(rvb, ab);
927 k = rvb0->wds - 1;
928 if (denorm)
929 /* do nothing */;
930 else if (rvb->wds <= k
931 || hi0bits( rvb->x[k]) >
932 hi0bits(rvb0->x[k])) {
933 /* unlikely; can only have lost 1 high bit */
934 if (rve1 == emin) {
935 --rvbits;
936 denorm = 1;
937 }
938 else {
939 rvb = lshift(rvb, 1);
940 --rve;
941 --rve1;
942 L = finished = 0;
943 }
944 }
945 }
946 else {
947 rvb = sum(rvb, ab);
948 k = rvb->wds - 1;
949 if (k >= rvb0->wds
950 || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
951 if (denorm) {
952 if (++rvbits == nbits)
953 denorm = 0;
954 }
955 else {
956 rshift(rvb, 1);
957 rve++;
958 rve1++;
959 L = 0;
960 }
961 }
962 }
963 Bfree(ab);
964 Bfree(rvb0);
965 if (finished)
966 break;
967
968 z = rve + rvbits;
969 if (y == z && L) {
970 /* Can we stop now? */
971 tol = dval(&adj) * 5e-16; /* > max rel error */
972 dval(&adj) = adj0 - .5;
973 if (dval(&adj) < -tol) {
974 if (adj0 > tol) {
975 irv |= inex;
976 break;
977 }
978 }
979 else if (dval(&adj) > tol && adj0 < 1. - tol) {
980 irv |= inex;
981 break;
982 }
983 }
984 bb0 = denorm ? 0 : trailz(rvb);
985 Bfree(bb);
986 Bfree(bd);
987 Bfree(bs);
988 Bfree(delta);
989 }
990 if (!denorm && (j = nbits - rvbits)) {
991 if (j > 0)
992 rvb = lshift(rvb, j);
993 else
994 rshift(rvb, -j);
995 rve -= j;
996 }
997 *exp = rve;
998 Bfree(bb);
999 Bfree(bd);
1000 Bfree(bs);
1001 Bfree(bd0);
1002 Bfree(delta);
1003 if (rve > fpi->emax) {
1004 switch(fpi->rounding & 3) {
1005 case FPI_Round_near:
1006 goto huge;
1007 case FPI_Round_up:
1008 if (!sign)
1009 goto huge;
1010 break;
1011 case FPI_Round_down:
1012 if (sign)
1013 goto huge;
1014 }
1015 /* Round to largest representable magnitude */
1016 Bfree(rvb);
1017 rvb = 0;
1018 irv = STRTOG_Normal | STRTOG_Inexlo;
1019 *exp = fpi->emax;
1020 b = bits;
1021 be = b + ((fpi->nbits + 31) >> 5);
1022 while(b < be)
1023 *b++ = -1;
1024 if ((j = fpi->nbits & 0x1f))
1025 *--be >>= (32 - j);
1026 goto ret;
1027 huge:
1028 rvb->wds = 0;
1029 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
1030 #ifndef NO_ERRNO
1031 errno = ERANGE;
1032 #endif
1033 infnanexp:
1034 *exp = fpi->emax + 1;
1035 }
1036 ret:
1037 if (denorm) {
1038 if (sudden_underflow) {
1039 rvb->wds = 0;
1040 irv = STRTOG_Underflow | STRTOG_Inexlo;
1041 #ifndef NO_ERRNO
1042 errno = ERANGE;
1043 #endif
1044 }
1045 else {
1046 irv = (irv & ~STRTOG_Retmask) |
1047 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
1048 if (irv & STRTOG_Inexact) {
1049 irv |= STRTOG_Underflow;
1050 #ifndef NO_ERRNO
1051 errno = ERANGE;
1052 #endif
1053 }
1054 }
1055 }
1056 if (se)
1057 *se = (char *)s;
1058 if (sign)
1059 irv |= STRTOG_Neg;
1060 if (rvb) {
1061 copybits(bits, nbits, rvb);
1062 Bfree(rvb);
1063 }
1064 return irv;
1065 }
1066