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