1 /* $NetBSD: strtodg.c,v 1.12 2013/04/19 10:41:53 joerg 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 (e < fpi->emin) {
254 k = fpi->emin - e;
255 e = fpi->emin;
256 if (k > nb || fpi->sudden_underflow) {
257 b->wds = inex = 0;
258 *irv = STRTOG_Underflow | STRTOG_Inexlo;
259 }
260 else {
261 k1 = k - 1;
262 if (k1 > 0 && !lostbits)
263 lostbits = any_on(b, k1);
264 if (!lostbits && !exact)
265 goto ret;
266 lostbits |=
267 carry = b->x[(unsigned int)k1>>kshift] &
268 (1 << (k1 & kmask));
269 rshift(b, k);
270 *irv = STRTOG_Denormal;
271 if (carry) {
272 b = increment(b);
273 inex = STRTOG_Inexhi | STRTOG_Underflow;
274 }
275 else if (lostbits)
276 inex = STRTOG_Inexlo | STRTOG_Underflow;
277 }
278 }
279 else if (e > fpi->emax) {
280 e = fpi->emax + 1;
281 *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
282 #ifndef NO_ERRNO
283 errno = ERANGE;
284 #endif
285 b->wds = inex = 0;
286 }
287 *expt = e;
288 copybits(bits, nb, b);
289 *irv |= inex;
290 rv = 1;
291 ret:
292 Bfree(b);
293 return rv;
294 }
295
296 #ifndef VAX
297 static int
298 #ifdef KR_headers
mantbits(d)299 mantbits(d) U *d;
300 #else
301 mantbits(U *d)
302 #endif
303 {
304 ULong L;
305 #ifdef VAX
306 L = word1(d) << 16 | word1(d) >> 16;
307 if (L)
308 #else
309 if ( (L = word1(d)) !=0)
310 #endif
311 return P - lo0bits(&L);
312 #ifdef VAX
313 L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;
314 #else
315 L = word0(d) | Exp_msk1;
316 #endif
317 return P - 32 - lo0bits(&L);
318 }
319 #endif /* !VAX */
320
321 int
strtodg(CONST char * s00,char ** se,CONST FPI * fpi,Long * expt,ULong * bits,locale_t loc)322 strtodg(CONST char *s00, char **se, CONST FPI *fpi, Long *expt, ULong *bits,
323 locale_t loc)
324 {
325 int abe, abits, asub;
326 #ifdef INFNAN_CHECK
327 int decpt;
328 #endif
329 int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, denorm;
330 int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
331 int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
332 int sudden_underflow = 0; /* pacify gcc */
333 CONST char *s, *s0, *s1;
334 double adj0, tol;
335 Long L;
336 U adj, rv;
337 ULong *b, *be, y, z;
338 Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
339 #ifdef USE_LOCALE /*{{*/
340 char *decimalpoint = localeconv_l(loc)->decimal_point;
341 size_t dplen = strlen(decimalpoint);
342 #endif /*USE_LOCALE}}*/
343
344 e2 = 0; /* XXX gcc */
345
346 irv = STRTOG_Zero;
347 denorm = sign = nz0 = nz = 0;
348 dval(&rv) = 0.;
349 rvb = 0;
350 nbits = fpi->nbits;
351 for(s = s00;;s++) switch(*s) {
352 case '-':
353 sign = 1;
354 /* FALLTHROUGH */
355 case '+':
356 if (*++s)
357 goto break2;
358 /* FALLTHROUGH */
359 case 0:
360 sign = 0;
361 irv = STRTOG_NoNumber;
362 s = s00;
363 goto ret;
364 case '\t':
365 case '\n':
366 case '\v':
367 case '\f':
368 case '\r':
369 case ' ':
370 continue;
371 default:
372 goto break2;
373 }
374 break2:
375 if (*s == '0') {
376 #ifndef NO_HEX_FP
377 switch(s[1]) {
378 case 'x':
379 case 'X':
380 irv = gethex(&s, fpi, expt, &rvb, sign, loc);
381 if (irv == STRTOG_NoNumber) {
382 s = s00;
383 sign = 0;
384 }
385 goto ret;
386 }
387 #endif
388 nz0 = 1;
389 while(*++s == '0') ;
390 if (!*s)
391 goto ret;
392 }
393 sudden_underflow = fpi->sudden_underflow;
394 s0 = s;
395 y = z = 0;
396 #ifdef INFNAN_CHECK
397 decpt = 0;
398 #endif
399 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
400 if (nd < 9)
401 y = 10*y + c - '0';
402 else if (nd < 16)
403 z = 10*z + c - '0';
404 nd0 = nd;
405 #ifdef USE_LOCALE
406 if (c == *decimalpoint) {
407 for(i = 1; decimalpoint[i]; ++i)
408 if (s[i] != decimalpoint[i])
409 goto dig_done;
410 s += i;
411 c = *s;
412 #else
413 if (c == '.') {
414 c = *++s;
415 #endif
416 #ifdef INFNAN_CHECK
417 decpt = 1;
418 #endif
419 if (!nd) {
420 for(; c == '0'; c = *++s)
421 nz++;
422 if (c > '0' && c <= '9') {
423 s0 = s;
424 nf += nz;
425 nz = 0;
426 goto have_dig;
427 }
428 goto dig_done;
429 }
430 for(; c >= '0' && c <= '9'; c = *++s) {
431 have_dig:
432 nz++;
433 if (c -= '0') {
434 nf += nz;
435 for(i = 1; i < nz; i++)
436 if (nd++ < 9)
437 y *= 10;
438 else if (nd <= DBL_DIG + 1)
439 z *= 10;
440 if (nd++ < 9)
441 y = 10*y + c;
442 else if (nd <= DBL_DIG + 1)
443 z = 10*z + c;
444 nz = 0;
445 }
446 }
447 }/*}*/
448 dig_done:
449 e = 0;
450 if (c == 'e' || c == 'E') {
451 if (!nd && !nz && !nz0) {
452 irv = STRTOG_NoNumber;
453 s = s00;
454 goto ret;
455 }
456 s00 = s;
457 esign = 0;
458 switch(c = *++s) {
459 case '-':
460 esign = 1;
461 /* FALLTHROUGH */
462 case '+':
463 c = *++s;
464 }
465 if (c >= '0' && c <= '9') {
466 while(c == '0')
467 c = *++s;
468 if (c > '0' && c <= '9') {
469 L = c - '0';
470 s1 = s;
471 while((c = *++s) >= '0' && c <= '9')
472 L = 10*L + c - '0';
473 if (s - s1 > 8 || L > 19999)
474 /* Avoid confusion from exponents
475 * so large that e might overflow.
476 */
477 e = 19999; /* safe for 16 bit ints */
478 else
479 e = (int)L;
480 if (esign)
481 e = -e;
482 }
483 else
484 e = 0;
485 }
486 else
487 s = s00;
488 }
489 if (!nd) {
490 if (!nz && !nz0) {
491 #ifdef INFNAN_CHECK
492 /* Check for Nan and Infinity */
493 if (!decpt)
494 switch(c) {
495 case 'i':
496 case 'I':
497 if (match(&s,"nf")) {
498 --s;
499 if (!match(&s,"inity"))
500 ++s;
501 irv = STRTOG_Infinite;
502 goto infnanexp;
503 }
504 break;
505 case 'n':
506 case 'N':
507 if (match(&s, "an")) {
508 irv = STRTOG_NaN;
509 *expt = fpi->emax + 1;
510 #ifndef No_Hex_NaN
511 if (*s == '(') /*)*/
512 irv = hexnan(&s, fpi, bits);
513 #endif
514 goto infnanexp;
515 }
516 }
517 #endif /* INFNAN_CHECK */
518 irv = STRTOG_NoNumber;
519 s = s00;
520 }
521 goto ret;
522 }
523
524 irv = STRTOG_Normal;
525 e1 = e -= nf;
526 rd = 0;
527 switch(fpi->rounding & 3) {
528 case FPI_Round_up:
529 rd = 2 - sign;
530 break;
531 case FPI_Round_zero:
532 rd = 1;
533 break;
534 case FPI_Round_down:
535 rd = 1 + sign;
536 }
537
538 /* Now we have nd0 digits, starting at s0, followed by a
539 * decimal point, followed by nd-nd0 digits. The number we're
540 * after is the integer represented by those digits times
541 * 10**e */
542
543 if (!nd0)
544 nd0 = nd;
545 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
546 dval(&rv) = y;
547 if (k > 9)
548 dval(&rv) = tens[k - 9] * dval(&rv) + z;
549 bd0 = 0;
550 if (nbits <= P && nd <= DBL_DIG) {
551 if (!e) {
552 if (rvOK(&rv, fpi, expt, bits, 1, rd, &irv))
553 goto ret;
554 }
555 else if (e > 0) {
556 if (e <= Ten_pmax) {
557 #ifdef VAX
558 goto vax_ovfl_check;
559 #else
560 i = fivesbits[e] + mantbits(&rv) <= P;
561 /* rv = */ rounded_product(dval(&rv), tens[e]);
562 if (rvOK(&rv, fpi, expt, bits, i, rd, &irv))
563 goto ret;
564 e1 -= e;
565 goto rv_notOK;
566 #endif
567 }
568 i = DBL_DIG - nd;
569 if (e <= Ten_pmax + i) {
570 /* A fancier test would sometimes let us do
571 * this for larger i values.
572 */
573 e2 = e - i;
574 e1 -= i;
575 dval(&rv) *= tens[i];
576 #ifdef VAX
577 /* VAX exponent range is so narrow we must
578 * worry about overflow here...
579 */
580 vax_ovfl_check:
581 dval(&adj) = dval(&rv);
582 word0(&adj) -= P*Exp_msk1;
583 /* adj = */ rounded_product(dval(&adj), tens[e2]);
584 if ((word0(&adj) & Exp_mask)
585 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
586 goto rv_notOK;
587 word0(&adj) += P*Exp_msk1;
588 dval(&rv) = dval(&adj);
589 #else
590 /* rv = */ rounded_product(dval(&rv), tens[e2]);
591 #endif
592 if (rvOK(&rv, fpi, expt, bits, 0, rd, &irv))
593 goto ret;
594 e1 -= e2;
595 }
596 }
597 #ifndef Inaccurate_Divide
598 else if (e >= -Ten_pmax) {
599 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
600 if (rvOK(&rv, fpi, expt, bits, 0, rd, &irv))
601 goto ret;
602 e1 -= e;
603 }
604 #endif
605 }
606 rv_notOK:
607 e1 += nd - k;
608
609 /* Get starting approximation = rv * 10**e1 */
610
611 e2 = 0;
612 if (e1 > 0) {
613 if ( (i = e1 & 15) !=0)
614 dval(&rv) *= tens[i];
615 if (e1 &= ~15) {
616 e1 = (unsigned int)e1 >> 4;
617 while(e1 >= (1 << (n_bigtens-1))) {
618 e2 += ((word0(&rv) & Exp_mask)
619 >> Exp_shift1) - Bias;
620 word0(&rv) &= ~Exp_mask;
621 word0(&rv) |= Bias << Exp_shift1;
622 dval(&rv) *= bigtens[n_bigtens-1];
623 e1 -= 1 << (n_bigtens-1);
624 }
625 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
626 word0(&rv) &= ~Exp_mask;
627 word0(&rv) |= Bias << Exp_shift1;
628 for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
629 if (e1 & 1)
630 dval(&rv) *= bigtens[j];
631 }
632 }
633 else if (e1 < 0) {
634 e1 = -e1;
635 if ( (i = e1 & 15) !=0)
636 dval(&rv) /= tens[i];
637 if (e1 &= ~15) {
638 e1 = (unsigned int)e1 >> 4;
639 while(e1 >= (1 << (n_bigtens-1))) {
640 e2 += ((word0(&rv) & Exp_mask)
641 >> Exp_shift1) - Bias;
642 word0(&rv) &= ~Exp_mask;
643 word0(&rv) |= Bias << Exp_shift1;
644 dval(&rv) *= tinytens[n_bigtens-1];
645 e1 -= 1 << (n_bigtens-1);
646 }
647 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
648 word0(&rv) &= ~Exp_mask;
649 word0(&rv) |= Bias << Exp_shift1;
650 for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
651 if (e1 & 1)
652 dval(&rv) *= tinytens[j];
653 }
654 }
655 #ifdef IBM
656 /* e2 is a correction to the (base 2) exponent of the return
657 * value, reflecting adjustments above to avoid overflow in the
658 * native arithmetic. For native IBM (base 16) arithmetic, we
659 * must multiply e2 by 4 to change from base 16 to 2.
660 */
661 e2 <<= 2;
662 #endif
663 rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */
664 if (rvb == NULL)
665 return STRTOG_NoMemory;
666 rve += e2;
667 if ((j = rvbits - nbits) > 0) {
668 rshift(rvb, j);
669 rvbits = nbits;
670 rve += j;
671 }
672 bb0 = 0; /* trailing zero bits in rvb */
673 e2 = rve + rvbits - nbits;
674 if (e2 > fpi->emax + 1)
675 goto huge;
676 rve1 = rve + rvbits - nbits;
677 if (e2 < (emin = fpi->emin)) {
678 denorm = 1;
679 j = rve - emin;
680 if (j > 0) {
681 rvb = lshift(rvb, j);
682 rvbits += j;
683 }
684 else if (j < 0) {
685 rvbits += j;
686 if (rvbits <= 0) {
687 if (rvbits < -1) {
688 ufl:
689 rvb->wds = 0;
690 rvb->x[0] = 0;
691 *expt = emin;
692 irv = STRTOG_Underflow | STRTOG_Inexlo;
693 goto ret;
694 }
695 rvb->x[0] = rvb->wds = rvbits = 1;
696 }
697 else
698 rshift(rvb, -j);
699 }
700 rve = rve1 = emin;
701 if (sudden_underflow && e2 + 1 < emin)
702 goto ufl;
703 }
704
705 /* Now the hard part -- adjusting rv to the correct value.*/
706
707 /* Put digits into bd: true value = bd * 10^e */
708
709 bd0 = s2b(s0, nd0, nd, y, dplen);
710
711 for(;;) {
712 bd = Balloc(bd0->k);
713 if (bd == NULL)
714 return STRTOG_NoMemory;
715 Bcopy(bd, bd0);
716 bb = Balloc(rvb->k);
717 if (bb == NULL)
718 return STRTOG_NoMemory;
719 Bcopy(bb, rvb);
720 bbbits = rvbits - bb0;
721 bbe = rve + bb0;
722 bs = i2b(1);
723 if (bs == NULL)
724 return STRTOG_NoMemory;
725
726 if (e >= 0) {
727 bb2 = bb5 = 0;
728 bd2 = bd5 = e;
729 }
730 else {
731 bb2 = bb5 = -e;
732 bd2 = bd5 = 0;
733 }
734 if (bbe >= 0)
735 bb2 += bbe;
736 else
737 bd2 -= bbe;
738 bs2 = bb2;
739 j = nbits + 1 - bbbits;
740 i = bbe + bbbits - nbits;
741 if (i < emin) /* denormal */
742 j += i - emin;
743 bb2 += j;
744 bd2 += j;
745 i = bb2 < bd2 ? bb2 : bd2;
746 if (i > bs2)
747 i = bs2;
748 if (i > 0) {
749 bb2 -= i;
750 bd2 -= i;
751 bs2 -= i;
752 }
753 if (bb5 > 0) {
754 bs = pow5mult(bs, bb5);
755 if (bs == NULL)
756 return STRTOG_NoMemory;
757 bb1 = mult(bs, bb);
758 if (bb1 == NULL)
759 return STRTOG_NoMemory;
760 Bfree(bb);
761 bb = bb1;
762 }
763 bb2 -= bb0;
764 if (bb2 > 0) {
765 bb = lshift(bb, bb2);
766 if (bb == NULL)
767 return STRTOG_NoMemory;
768 }
769 else if (bb2 < 0)
770 rshift(bb, -bb2);
771 if (bd5 > 0) {
772 bd = pow5mult(bd, bd5);
773 if (bd == NULL)
774 return STRTOG_NoMemory;
775 }
776 if (bd2 > 0) {
777 bd = lshift(bd, bd2);
778 if (bd == NULL)
779 return STRTOG_NoMemory;
780 }
781 if (bs2 > 0) {
782 bs = lshift(bs, bs2);
783 if (bs == NULL)
784 return STRTOG_NoMemory;
785 }
786 asub = 1;
787 inex = STRTOG_Inexhi;
788 delta = diff(bb, bd);
789 if (delta == NULL)
790 return STRTOG_NoMemory;
791 if (delta->wds <= 1 && !delta->x[0])
792 break;
793 dsign = delta->sign;
794 delta->sign = finished = 0;
795 L = 0;
796 i = cmp(delta, bs);
797 if (rd && i <= 0) {
798 irv = STRTOG_Normal;
799 if ( (finished = dsign ^ (rd&1)) !=0) {
800 if (dsign != 0) {
801 irv |= STRTOG_Inexhi;
802 goto adj1;
803 }
804 irv |= STRTOG_Inexlo;
805 if (rve1 == emin)
806 goto adj1;
807 for(i = 0, j = nbits; j >= ULbits;
808 i++, j -= ULbits) {
809 if (rvb->x[i] & ALL_ON)
810 goto adj1;
811 }
812 if (j > 1 && lo0bits(rvb->x + i) < j - 1)
813 goto adj1;
814 rve = rve1 - 1;
815 rvb = set_ones(rvb, rvbits = nbits);
816 if (rvb == NULL)
817 return STRTOG_NoMemory;
818 break;
819 }
820 irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
821 break;
822 }
823 if (i < 0) {
824 /* Error is less than half an ulp -- check for
825 * special case of mantissa a power of two.
826 */
827 irv = dsign
828 ? STRTOG_Normal | STRTOG_Inexlo
829 : STRTOG_Normal | STRTOG_Inexhi;
830 if (dsign || bbbits > 1 || denorm || rve1 == emin)
831 break;
832 delta = lshift(delta,1);
833 if (delta == NULL)
834 return STRTOG_NoMemory;
835 if (cmp(delta, bs) > 0) {
836 irv = STRTOG_Normal | STRTOG_Inexlo;
837 goto drop_down;
838 }
839 break;
840 }
841 if (i == 0) {
842 /* exactly half-way between */
843 if (dsign) {
844 if (denorm && all_on(rvb, rvbits)) {
845 /*boundary case -- increment exponent*/
846 rvb->wds = 1;
847 rvb->x[0] = 1;
848 rve = emin + nbits - (rvbits = 1);
849 irv = STRTOG_Normal | STRTOG_Inexhi;
850 denorm = 0;
851 break;
852 }
853 irv = STRTOG_Normal | STRTOG_Inexlo;
854 }
855 else if (bbbits == 1) {
856 irv = STRTOG_Normal;
857 drop_down:
858 /* boundary case -- decrement exponent */
859 if (rve1 == emin) {
860 irv = STRTOG_Normal | STRTOG_Inexhi;
861 if (rvb->wds == 1 && rvb->x[0] == 1)
862 sudden_underflow = 1;
863 break;
864 }
865 rve -= nbits;
866 rvb = set_ones(rvb, rvbits = nbits);
867 if (rvb == NULL)
868 return STRTOG_NoMemory;
869 break;
870 }
871 else
872 irv = STRTOG_Normal | STRTOG_Inexhi;
873 if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
874 break;
875 if (dsign) {
876 rvb = increment(rvb);
877 if (rvb == NULL)
878 return STRTOG_NoMemory;
879 j = kmask & (ULbits - (rvbits & kmask));
880 if (hi0bits(rvb->x[rvb->wds - 1]) != j)
881 rvbits++;
882 irv = STRTOG_Normal | STRTOG_Inexhi;
883 }
884 else {
885 if (bbbits == 1)
886 goto undfl;
887 decrement(rvb);
888 irv = STRTOG_Normal | STRTOG_Inexlo;
889 }
890 break;
891 }
892 if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
893 adj1:
894 inex = STRTOG_Inexlo;
895 if (dsign) {
896 asub = 0;
897 inex = STRTOG_Inexhi;
898 }
899 else if (denorm && bbbits <= 1) {
900 undfl:
901 rvb->wds = 0;
902 rve = emin;
903 irv = STRTOG_Underflow | STRTOG_Inexlo;
904 break;
905 }
906 adj0 = dval(&adj) = 1.;
907 }
908 else {
909 adj0 = dval(&adj) *= 0.5;
910 if (dsign) {
911 asub = 0;
912 inex = STRTOG_Inexlo;
913 }
914 if (dval(&adj) < 2147483647.) {
915 L = adj0;
916 adj0 -= L;
917 switch(rd) {
918 case 0:
919 if (adj0 >= .5)
920 goto inc_L;
921 break;
922 case 1:
923 if (asub && adj0 > 0.)
924 goto inc_L;
925 break;
926 case 2:
927 if (!asub && adj0 > 0.) {
928 inc_L:
929 L++;
930 inex = STRTOG_Inexact - inex;
931 }
932 }
933 dval(&adj) = L;
934 }
935 }
936 y = rve + rvbits;
937
938 /* adj *= ulp(dval(&rv)); */
939 /* if (asub) rv -= adj; else rv += adj; */
940
941 if (!denorm && rvbits < nbits) {
942 rvb = lshift(rvb, j = nbits - rvbits);
943 if (rvb == NULL)
944 return STRTOG_NoMemory;
945 rve -= j;
946 rvbits = nbits;
947 }
948 ab = d2b(dval(&adj), &abe, &abits);
949 if (ab == NULL)
950 return STRTOG_NoMemory;
951 if (abe < 0)
952 rshift(ab, -abe);
953 else if (abe > 0)
954 ab = lshift(ab, abe);
955 rvb0 = rvb;
956 if (asub) {
957 /* rv -= adj; */
958 j = hi0bits(rvb->x[rvb->wds-1]);
959 rvb = diff(rvb, ab);
960 if (rvb == NULL)
961 return STRTOG_NoMemory;
962 k = rvb0->wds - 1;
963 if (denorm)
964 /* do nothing */;
965 else if (rvb->wds <= k
966 || hi0bits( rvb->x[k]) >
967 hi0bits(rvb0->x[k])) {
968 /* unlikely; can only have lost 1 high bit */
969 if (rve1 == emin) {
970 --rvbits;
971 denorm = 1;
972 }
973 else {
974 rvb = lshift(rvb, 1);
975 if (rvb == NULL)
976 return STRTOG_NoMemory;
977 --rve;
978 --rve1;
979 L = finished = 0;
980 }
981 }
982 }
983 else {
984 rvb = sum(rvb, ab);
985 if (rvb == NULL)
986 return STRTOG_NoMemory;
987 k = rvb->wds - 1;
988 if (k >= rvb0->wds
989 || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
990 if (denorm) {
991 if (++rvbits == nbits)
992 denorm = 0;
993 }
994 else {
995 rshift(rvb, 1);
996 rve++;
997 rve1++;
998 L = 0;
999 }
1000 }
1001 }
1002 Bfree(ab);
1003 Bfree(rvb0);
1004 if (finished)
1005 break;
1006
1007 z = rve + rvbits;
1008 if (y == z && L) {
1009 /* Can we stop now? */
1010 tol = dval(&adj) * 5e-16; /* > max rel error */
1011 dval(&adj) = adj0 - .5;
1012 if (dval(&adj) < -tol) {
1013 if (adj0 > tol) {
1014 irv |= inex;
1015 break;
1016 }
1017 }
1018 else if (dval(&adj) > tol && adj0 < 1. - tol) {
1019 irv |= inex;
1020 break;
1021 }
1022 }
1023 bb0 = denorm ? 0 : trailz(rvb);
1024 Bfree(bb);
1025 Bfree(bd);
1026 Bfree(bs);
1027 Bfree(delta);
1028 }
1029 if (!denorm && (j = nbits - rvbits)) {
1030 if (j > 0)
1031 rvb = lshift(rvb, j);
1032 else
1033 rshift(rvb, -j);
1034 rve -= j;
1035 }
1036 *expt = rve;
1037 Bfree(bb);
1038 Bfree(bd);
1039 Bfree(bs);
1040 Bfree(bd0);
1041 Bfree(delta);
1042 if (rve > fpi->emax) {
1043 switch(fpi->rounding & 3) {
1044 case FPI_Round_near:
1045 goto huge;
1046 case FPI_Round_up:
1047 if (!sign)
1048 goto huge;
1049 break;
1050 case FPI_Round_down:
1051 if (sign)
1052 goto huge;
1053 }
1054 /* Round to largest representable magnitude */
1055 Bfree(rvb);
1056 rvb = 0;
1057 irv = STRTOG_Normal | STRTOG_Inexlo;
1058 *expt = fpi->emax;
1059 b = bits;
1060 be = b + ((unsigned int)(fpi->nbits + 31) >> 5);
1061 while(b < be)
1062 *b++ = (unsigned int)-1;
1063 if ((j = fpi->nbits & 0x1f) != 0)
1064 *--be >>= (32 - j);
1065 goto ret;
1066 huge:
1067 rvb->wds = 0;
1068 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
1069 #ifndef NO_ERRNO
1070 errno = ERANGE;
1071 #endif
1072 #ifdef INFNAN_CHECK
1073 infnanexp:
1074 #endif
1075 *expt = fpi->emax + 1;
1076 }
1077 ret:
1078 if (denorm) {
1079 if (sudden_underflow) {
1080 rvb->wds = 0;
1081 irv = STRTOG_Underflow | STRTOG_Inexlo;
1082 #ifndef NO_ERRNO
1083 errno = ERANGE;
1084 #endif
1085 }
1086 else {
1087 irv = (irv & ~STRTOG_Retmask) |
1088 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
1089 if (irv & STRTOG_Inexact) {
1090 irv |= STRTOG_Underflow;
1091 #ifndef NO_ERRNO
1092 errno = ERANGE;
1093 #endif
1094 }
1095 }
1096 }
1097 if (se)
1098 *se = __UNCONST(s);
1099 if (sign)
1100 irv |= STRTOG_Neg;
1101 if (rvb) {
1102 copybits(bits, nbits, rvb);
1103 Bfree(rvb);
1104 }
1105 return irv;
1106 }
1107