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