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