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