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