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