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