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