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