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