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 "gdtoaimp.h"
33
34#ifdef USE_LOCALE
35#include "locale.h"
36#endif
37
38static const int
39fivesbits[] = { 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
44Bigint *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
84void 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
110static 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
124Bigint *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
147static 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
260static 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
269int __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, j, k;
274 int 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 /* fallthrough */
313 case '+':
314 if (*++s)
315 goto break2;
316 /* fallthrough */
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 /* fallthrough */
415 case '+':
416 c = *++s;
417 /* fallthrough */
418 }
419 if (c >= '0' && c <= '9') {
420 while(c == '0')
421 c = *++s;
422 if (c > '0' && c <= '9') {
423 L = c - '0';
424 s1 = s;
425 while((c = *++s) >= '0' && c <= '9')
426 L = 10*L + c - '0';
427 if (s - s1 > 8 || L > 19999)
428 /* Avoid confusion from exponents
429 * so large that e might overflow.
430 */
431 e = 19999; /* safe for 16 bit ints */
432 else
433 e = (int)L;
434 if (esign)
435 e = -e;
436 }
437 else
438 e = 0;
439 }
440 else
441 s = s00;
442 }
443 if (!nd) {
444 if (!nz && !nz0) {
445#ifdef INFNAN_CHECK
446 /* Check for Nan and Infinity */
447 if (!decpt)
448 switch(c) {
449 case 'i':
450 case 'I':
451 if (match(&s,"nf")) {
452 --s;
453 if (!match(&s,"inity"))
454 ++s;
455 irv = STRTOG_Infinite;
456 goto infnanexp;
457 }
458 break;
459 case 'n':
460 case 'N':
461 if (match(&s, "an")) {
462 irv = STRTOG_NaN;
463 *expo = fpi->emax + 1;
464#ifndef No_Hex_NaN
465 if (*s == '(') /*)*/
466 irv = hexnan(&s, fpi, bits);
467#endif
468 goto infnanexp;
469 }
470 }
471#endif /* INFNAN_CHECK */
472 irv = STRTOG_NoNumber;
473 s = s00;
474 }
475 goto ret;
476 }
477
478 irv = STRTOG_Normal;
479 e1 = e -= nf;
480 rd = 0;
481 switch(fpi->rounding & 3) {
482 case FPI_Round_up:
483 rd = 2 - sign;
484 break;
485 case FPI_Round_zero:
486 rd = 1;
487 break;
488 case FPI_Round_down:
489 rd = 1 + sign;
490 }
491
492 /* Now we have nd0 digits, starting at s0, followed by a
493 * decimal point, followed by nd-nd0 digits. The number we're
494 * after is the integer represented by those digits times
495 * 10**e */
496
497 if (!nd0)
498 nd0 = nd;
499 k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
500 dval(&rv) = y;
501 if (k > 9)
502 dval(&rv) = tens[k - 9] * dval(&rv) + z;
503 bd0 = 0;
504 if (nbits <= P && nd <= DBL_DIG) {
505 if (!e) {
506 if (rvOK(&rv, fpi, expo, bits, 1, rd, &irv))
507 goto ret;
508 }
509 else if (e > 0) {
510 if (e <= Ten_pmax) {
511 i = fivesbits[e] + mantbits(&rv) <= P;
512 /* rv = */ rounded_product(dval(&rv), tens[e]);
513 if (rvOK(&rv, fpi, expo, bits, i, rd, &irv))
514 goto ret;
515 e1 -= e;
516 goto rv_notOK;
517 }
518 i = DBL_DIG - nd;
519 if (e <= Ten_pmax + i) {
520 /* A fancier test would sometimes let us do
521 * this for larger i values.
522 */
523 e2 = e - i;
524 e1 -= i;
525 dval(&rv) *= tens[i];
526 /* rv = */ rounded_product(dval(&rv), tens[e2]);
527 if (rvOK(&rv, fpi, expo, bits, 0, rd, &irv))
528 goto ret;
529 e1 -= e2;
530 }
531 }
532#ifndef Inaccurate_Divide
533 else if (e >= -Ten_pmax) {
534 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
535 if (rvOK(&rv, fpi, expo, bits, 0, rd, &irv))
536 goto ret;
537 e1 -= e;
538 }
539#endif
540 }
541 rv_notOK:
542 e1 += nd - k;
543
544 /* Get starting approximation = rv * 10**e1 */
545
546 e2 = 0;
547 if (e1 > 0) {
548 if ( (i = e1 & 15) !=0)
549 dval(&rv) *= tens[i];
550 if (e1 &= ~15) {
551 e1 >>= 4;
552 while(e1 >= (1 << (n_bigtens-1))) {
553 e2 += ((word0(&rv) & Exp_mask)
554 >> Exp_shift1) - Bias;
555 word0(&rv) &= ~Exp_mask;
556 word0(&rv) |= Bias << Exp_shift1;
557 dval(&rv) *= bigtens[n_bigtens-1];
558 e1 -= 1 << (n_bigtens-1);
559 }
560 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
561 word0(&rv) &= ~Exp_mask;
562 word0(&rv) |= Bias << Exp_shift1;
563 for(j = 0; e1 > 0; j++, e1 >>= 1)
564 if (e1 & 1)
565 dval(&rv) *= bigtens[j];
566 }
567 }
568 else if (e1 < 0) {
569 e1 = -e1;
570 if ( (i = e1 & 15) !=0)
571 dval(&rv) /= tens[i];
572 if (e1 &= ~15) {
573 e1 >>= 4;
574 while(e1 >= (1 << (n_bigtens-1))) {
575 e2 += ((word0(&rv) & Exp_mask)
576 >> Exp_shift1) - Bias;
577 word0(&rv) &= ~Exp_mask;
578 word0(&rv) |= Bias << Exp_shift1;
579 dval(&rv) *= tinytens[n_bigtens-1];
580 e1 -= 1 << (n_bigtens-1);
581 }
582 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
583 word0(&rv) &= ~Exp_mask;
584 word0(&rv) |= Bias << Exp_shift1;
585 for(j = 0; e1 > 0; j++, e1 >>= 1)
586 if (e1 & 1)
587 dval(&rv) *= tinytens[j];
588 }
589 }
590 rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */
591 rve += e2;
592 if ((j = rvbits - nbits) > 0) {
593 rshift(rvb, j);
594 rvbits = nbits;
595 rve += j;
596 }
597 bb0 = 0; /* trailing zero bits in rvb */
598 e2 = rve + rvbits - nbits;
599 if (e2 > fpi->emax + 1)
600 goto huge;
601 rve1 = rve + rvbits - nbits;
602 if (e2 < (emin = fpi->emin)) {
603 denorm = 1;
604 j = rve - emin;
605 if (j > 0) {
606 rvb = lshift(rvb, j);
607 rvbits += j;
608 }
609 else if (j < 0) {
610 rvbits += j;
611 if (rvbits <= 0) {
612 if (rvbits < -1) {
613 ufl:
614 rvb->wds = 0;
615 rvb->x[0] = 0;
616 *expo = emin;
617 irv = STRTOG_Underflow | STRTOG_Inexlo;
618 goto ret;
619 }
620 rvb->x[0] = rvb->wds = rvbits = 1;
621 }
622 else
623 rshift(rvb, -j);
624 }
625 rve = rve1 = emin;
626 if (sudden_underflow && e2 + 1 < emin)
627 goto ufl;
628 }
629
630 /* Now the hard part -- adjusting rv to the correct value.*/
631
632 /* Put digits into bd: true value = bd * 10^e */
633
634 bd0 = s2b(s0, nd0, nd, y, dplen);
635
636 for(;;) {
637 bd = Balloc(bd0->k);
638 Bcopy(bd, bd0);
639 bb = Balloc(rvb->k);
640 Bcopy(bb, rvb);
641 bbbits = rvbits - bb0;
642 bbe = rve + bb0;
643 bs = i2b(1);
644
645 if (e >= 0) {
646 bb2 = bb5 = 0;
647 bd2 = bd5 = e;
648 }
649 else {
650 bb2 = bb5 = -e;
651 bd2 = bd5 = 0;
652 }
653 if (bbe >= 0)
654 bb2 += bbe;
655 else
656 bd2 -= bbe;
657 bs2 = bb2;
658 j = nbits + 1 - bbbits;
659 i = bbe + bbbits - nbits;
660 if (i < emin) /* denormal */
661 j += i - emin;
662 bb2 += j;
663 bd2 += j;
664 i = bb2 < bd2 ? bb2 : bd2;
665 if (i > bs2)
666 i = bs2;
667 if (i > 0) {
668 bb2 -= i;
669 bd2 -= i;
670 bs2 -= i;
671 }
672 if (bb5 > 0) {
673 bs = pow5mult(bs, bb5);
674 bb1 = mult(bs, bb);
675 Bfree(bb);
676 bb = bb1;
677 }
678 bb2 -= bb0;
679 if (bb2 > 0)
680 bb = lshift(bb, bb2);
681 else if (bb2 < 0)
682 rshift(bb, -bb2);
683 if (bd5 > 0)
684 bd = pow5mult(bd, bd5);
685 if (bd2 > 0)
686 bd = lshift(bd, bd2);
687 if (bs2 > 0)
688 bs = lshift(bs, bs2);
689 asub = 1;
690 inex = STRTOG_Inexhi;
691 delta = diff(bb, bd);
692 if (delta->wds <= 1 && !delta->x[0])
693 break;
694 dsign = delta->sign;
695 delta->sign = finished = 0;
696 L = 0;
697 i = cmp(delta, bs);
698 if (rd && i <= 0) {
699 irv = STRTOG_Normal;
700 if ( (finished = dsign ^ (rd&1)) !=0) {
701 if (dsign != 0) {
702 irv |= STRTOG_Inexhi;
703 goto adj1;
704 }
705 irv |= STRTOG_Inexlo;
706 if (rve1 == emin)
707 goto adj1;
708 for(i = 0, j = nbits; j >= ULbits;
709 i++, j -= ULbits) {
710 if (rvb->x[i] & ALL_ON)
711 goto adj1;
712 }
713 if (j > 1 && lo0bits(rvb->x + i) < j - 1)
714 goto adj1;
715 rve = rve1 - 1;
716 rvb = set_ones(rvb, rvbits = nbits);
717 break;
718 }
719 irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
720 break;
721 }
722 if (i < 0) {
723 /* Error is less than half an ulp -- check for
724 * special case of mantissa a power of two.
725 */
726 irv = dsign
727 ? STRTOG_Normal | STRTOG_Inexlo
728 : STRTOG_Normal | STRTOG_Inexhi;
729 if (dsign || bbbits > 1 || denorm || rve1 == emin)
730 break;
731 delta = lshift(delta,1);
732 if (cmp(delta, bs) > 0) {
733 irv = STRTOG_Normal | STRTOG_Inexlo;
734 goto drop_down;
735 }
736 break;
737 }
738 if (i == 0) {
739 /* exactly half-way between */
740 if (dsign) {
741 if (denorm && all_on(rvb, rvbits)) {
742 /*boundary case -- increment exponent*/
743 rvb->wds = 1;
744 rvb->x[0] = 1;
745 rve = emin + nbits - (rvbits = 1);
746 irv = STRTOG_Normal | STRTOG_Inexhi;
747 denorm = 0;
748 break;
749 }
750 irv = STRTOG_Normal | STRTOG_Inexlo;
751 }
752 else if (bbbits == 1) {
753 irv = STRTOG_Normal;
754 drop_down:
755 /* boundary case -- decrement exponent */
756 if (rve1 == emin) {
757 irv = STRTOG_Normal | STRTOG_Inexhi;
758 if (rvb->wds == 1 && rvb->x[0] == 1)
759 sudden_underflow = 1;
760 break;
761 }
762 rve -= nbits;
763 rvb = set_ones(rvb, rvbits = nbits);
764 break;
765 }
766 else
767 irv = STRTOG_Normal | STRTOG_Inexhi;
768 if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
769 break;
770 if (dsign) {
771 rvb = increment(rvb);
772 j = kmask & (ULbits - (rvbits & kmask));
773 if (hi0bits(rvb->x[rvb->wds - 1]) != j)
774 rvbits++;
775 irv = STRTOG_Normal | STRTOG_Inexhi;
776 }
777 else {
778 if (bbbits == 1)
779 goto undfl;
780 decrement(rvb);
781 irv = STRTOG_Normal | STRTOG_Inexlo;
782 }
783 break;
784 }
785 if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
786 adj1:
787 inex = STRTOG_Inexlo;
788 if (dsign) {
789 asub = 0;
790 inex = STRTOG_Inexhi;
791 }
792 else if (denorm && bbbits <= 1) {
793 undfl:
794 rvb->wds = 0;
795 rve = emin;
796 irv = STRTOG_Underflow | STRTOG_Inexlo;
797 break;
798 }
799 adj0 = dval(&adj) = 1.;
800 }
801 else {
802 adj0 = dval(&adj) *= 0.5;
803 if (dsign) {
804 asub = 0;
805 inex = STRTOG_Inexlo;
806 }
807 if (dval(&adj) < 2147483647.) {
808 L = adj0;
809 adj0 -= L;
810 switch(rd) {
811 case 0:
812 if (adj0 >= .5)
813 goto inc_L;
814 break;
815 case 1:
816 if (asub && adj0 > 0.)
817 goto inc_L;
818 break;
819 case 2:
820 if (!asub && adj0 > 0.) {
821 inc_L:
822 L++;
823 inex = STRTOG_Inexact - inex;
824 }
825 }
826 dval(&adj) = L;
827 }
828 }
829 y = rve + rvbits;
830
831 /* adj *= ulp(&rv); */
832 /* if (asub) rv -= adj; else rv += adj; */
833
834 if (!denorm && rvbits < nbits) {
835 rvb = lshift(rvb, j = nbits - rvbits);
836 rve -= j;
837 rvbits = nbits;
838 }
839 ab = d2b(dval(&adj), &abe, &abits);
840 if (abe < 0)
841 rshift(ab, -abe);
842 else if (abe > 0)
843 ab = lshift(ab, abe);
844 rvb0 = rvb;
845 if (asub) {
846 /* rv -= adj; */
847 j = hi0bits(rvb->x[rvb->wds-1]);
848 rvb = diff(rvb, ab);
849 k = rvb0->wds - 1;
850 if (denorm)
851 /* do nothing */;
852 else if (rvb->wds <= k
853 || hi0bits( rvb->x[k]) >
854 hi0bits(rvb0->x[k])) {
855 /* unlikely; can only have lost 1 high bit */
856 if (rve1 == emin) {
857 --rvbits;
858 denorm = 1;
859 }
860 else {
861 rvb = lshift(rvb, 1);
862 --rve;
863 --rve1;
864 L = finished = 0;
865 }
866 }
867 }
868 else {
869 rvb = sum(rvb, ab);
870 k = rvb->wds - 1;
871 if (k >= rvb0->wds
872 || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
873 if (denorm) {
874 if (++rvbits == nbits)
875 denorm = 0;
876 }
877 else {
878 rshift(rvb, 1);
879 rve++;
880 rve1++;
881 L = 0;
882 }
883 }
884 }
885 Bfree(ab);
886 Bfree(rvb0);
887 if (finished)
888 break;
889
890 z = rve + rvbits;
891 if (y == z && L) {
892 /* Can we stop now? */
893 tol = dval(&adj) * 5e-16; /* > max rel error */
894 dval(&adj) = adj0 - .5;
895 if (dval(&adj) < -tol) {
896 if (adj0 > tol) {
897 irv |= inex;
898 break;
899 }
900 }
901 else if (dval(&adj) > tol && adj0 < 1. - tol) {
902 irv |= inex;
903 break;
904 }
905 }
906 bb0 = denorm ? 0 : trailz(rvb);
907 Bfree(bb);
908 Bfree(bd);
909 Bfree(bs);
910 Bfree(delta);
911 }
912 if (!denorm && (j = nbits - rvbits)) {
913 if (j > 0)
914 rvb = lshift(rvb, j);
915 else
916 rshift(rvb, -j);
917 rve -= j;
918 }
919 *expo = rve;
920 Bfree(bb);
921 Bfree(bd);
922 Bfree(bs);
923 Bfree(bd0);
924 Bfree(delta);
925 if (rve > fpi->emax) {
926huge:
927 Bfree(rvb);
928 rvb = 0;
929 SET_ERRNO(ERANGE);
930 switch(fpi->rounding & 3) {
931 case FPI_Round_up:
932 if (!sign)
933 goto ret_inf;
934 break;
935 case FPI_Round_down:
936 if (!sign)
937 break;
938 /* fallthrough */
939 case FPI_Round_near:
940 ret_inf:
941 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
942 k = nbits >> kshift;
943 if (nbits & kmask)
944 ++k;
945 memset(bits, 0, k*sizeof(ULong));
946 infnanexp:
947 *expo = fpi->emax + 1;
948 goto ret;
949 }
950 /* Round to largest representable magnitude */
951 irv = STRTOG_Normal | STRTOG_Inexlo;
952 *expo = fpi->emax;
953 b = bits;
954 be = b + ((fpi->nbits + 31) >> 5);
955 while(b < be)
956 *b++ = -1;
957 if ((j = fpi->nbits & 0x1f))
958 *--be >>= (32 - j);
959 }
960 ret:
961 if (denorm) {
962 if (sudden_underflow) {
963 rvb->wds = 0;
964 irv = STRTOG_Underflow | STRTOG_Inexlo;
965 SET_ERRNO(ERANGE);
966 }
967 else {
968 irv = (irv & ~STRTOG_Retmask) |
969 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
970 if (irv & STRTOG_Inexact) {
971 irv |= STRTOG_Underflow;
972 SET_ERRNO(ERANGE);
973 }
974 }
975 }
976 if (se)
977 *se = (char *)s;
978 if (sign)
979 irv |= STRTOG_Neg;
980 if (rvb) {
981 copybits(bits, nbits, rvb);
982 Bfree(rvb);
983 }
984 return irv;
985}