]> git.saurik.com Git - apple/libc.git/blob - gdtoa/FreeBSD/gdtoa-dtoa.c
Libc-1439.100.3.tar.gz
[apple/libc.git] / gdtoa / FreeBSD / gdtoa-dtoa.c
1 /****************************************************************
2
3 The author of this software is David M. Gay.
4
5 Copyright (C) 1998, 1999 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 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
35 *
36 * Inspired by "How to Print Floating-Point Numbers Accurately" by
37 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
38 *
39 * Modifications:
40 * 1. Rather than iterating, we use a simple numeric overestimate
41 * to determine k = floor(log10(d)). We scale relevant
42 * quantities using O(log2(k)) rather than O(k) multiplications.
43 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
44 * try to generate digits strictly left to right. Instead, we
45 * compute with fewer bits and propagate the carry if necessary
46 * when rounding the final digit up. This is often faster.
47 * 3. Under the assumption that input will be rounded nearest,
48 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
49 * That is, we allow equality in stopping tests when the
50 * round-nearest rule will give the same floating-point value
51 * as would satisfaction of the stopping test with strict
52 * inequality.
53 * 4. We remove common factors of powers of 2 from relevant
54 * quantities.
55 * 5. When converting floating-point integers less than 1e16,
56 * we use floating-point arithmetic rather than resorting
57 * to multiple-precision integers.
58 * 6. When asked to produce fewer than 15 digits, we first try
59 * to get by with floating-point arithmetic; we resort to
60 * multiple-precision integer arithmetic only if we cannot
61 * guarantee that the floating-point calculation has given
62 * the correctly rounded result. For k requested digits and
63 * "uniformly" distributed input, the probability is
64 * something like 10^(k-15) that we must resort to the Long
65 * calculation.
66 */
67
68 #ifdef Honor_FLT_ROUNDS
69 #undef Check_FLT_ROUNDS
70 #define Check_FLT_ROUNDS
71 #else
72 #define Rounding Flt_Rounds
73 #endif
74
75 char *
76 dtoa
77 #ifdef KR_headers
78 (d0, mode, ndigits, decpt, sign, rve)
79 double d0; int mode, ndigits, *decpt, *sign; char **rve;
80 #else
81 (double d0, int mode, int ndigits, int *decpt, int *sign, char **rve)
82 #endif
83 {
84 /* Arguments ndigits, decpt, sign are similar to those
85 of ecvt and fcvt; trailing zeros are suppressed from
86 the returned string. If not null, *rve is set to point
87 to the end of the return value. If d is +-Infinity or NaN,
88 then *decpt is set to 9999.
89
90 mode:
91 0 ==> shortest string that yields d when read in
92 and rounded to nearest.
93 1 ==> like 0, but with Steele & White stopping rule;
94 e.g. with IEEE P754 arithmetic , mode 0 gives
95 1e23 whereas mode 1 gives 9.999999999999999e22.
96 2 ==> max(1,ndigits) significant digits. This gives a
97 return value similar to that of ecvt, except
98 that trailing zeros are suppressed.
99 3 ==> through ndigits past the decimal point. This
100 gives a return value similar to that from fcvt,
101 except that trailing zeros are suppressed, and
102 ndigits can be negative.
103 4,5 ==> similar to 2 and 3, respectively, but (in
104 round-nearest mode) with the tests of mode 0 to
105 possibly return a shorter string that rounds to d.
106 With IEEE arithmetic and compilation with
107 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
108 as modes 2 and 3 when FLT_ROUNDS != 1.
109 6-9 ==> Debugging modes similar to mode - 4: don't try
110 fast floating-point estimate (if applicable).
111
112 Values of mode other than 0-9 are treated as mode 0.
113
114 Sufficient space is allocated to the return value
115 to hold the suppressed trailing zeros.
116 */
117
118 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
119 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
120 spec_case, try_quick;
121 Long L;
122 #ifndef Sudden_Underflow
123 int denorm;
124 ULong x;
125 #endif
126 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
127 U d, d2, eps;
128 double ds;
129 char *s, *s0;
130 #ifdef SET_INEXACT
131 int inexact, oldinexact;
132 #endif
133 #ifdef Honor_FLT_ROUNDS /*{*/
134 int Rounding;
135 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
136 Rounding = Flt_Rounds;
137 #else /*}{*/
138 Rounding = 1;
139 switch(fegetround()) {
140 case FE_TOWARDZERO: Rounding = 0; break;
141 case FE_UPWARD: Rounding = 2; break;
142 case FE_DOWNWARD: Rounding = 3;
143 }
144 #endif /*}}*/
145 #endif /*}*/
146
147 #ifndef MULTIPLE_THREADS
148 if (dtoa_result) {
149 freedtoa(dtoa_result);
150 dtoa_result = 0;
151 }
152 #endif
153 d.d = d0;
154 if (word0(&d) & Sign_bit) {
155 /* set sign for everything, including 0's and NaNs */
156 *sign = 1;
157 word0(&d) &= ~Sign_bit; /* clear sign bit */
158 }
159 else
160 *sign = 0;
161
162 #if defined(IEEE_Arith) + defined(VAX)
163 #ifdef IEEE_Arith
164 if ((word0(&d) & Exp_mask) == Exp_mask)
165 #else
166 if (word0(&d) == 0x8000)
167 #endif
168 {
169 /* Infinity or NaN */
170 *decpt = 9999;
171 #ifdef IEEE_Arith
172 if (!word1(&d) && !(word0(&d) & 0xfffff))
173 return nrv_alloc("Infinity", rve, 8);
174 #endif
175 return nrv_alloc("NaN", rve, 3);
176 }
177 #endif
178 #ifdef IBM
179 dval(&d) += 0; /* normalize */
180 #endif
181 if (!dval(&d)) {
182 *decpt = 1;
183 return nrv_alloc("0", rve, 1);
184 }
185
186 #ifdef SET_INEXACT
187 try_quick = oldinexact = get_inexact();
188 inexact = 1;
189 #endif
190 #ifdef Honor_FLT_ROUNDS
191 if (Rounding >= 2) {
192 if (*sign)
193 Rounding = Rounding == 2 ? 0 : 2;
194 else
195 if (Rounding != 2)
196 Rounding = 0;
197 }
198 #endif
199
200 b = d2b(dval(&d), &be, &bbits);
201 #ifdef Sudden_Underflow
202 i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
203 #else
204 if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
205 #endif
206 dval(&d2) = dval(&d);
207 word0(&d2) &= Frac_mask1;
208 word0(&d2) |= Exp_11;
209 #ifdef IBM
210 if (( j = 11 - hi0bits(word0(&d2) & Frac_mask) )!=0)
211 dval(&d2) /= 1 << j;
212 #endif
213
214 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
215 * log10(x) = log(x) / log(10)
216 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
217 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
218 *
219 * This suggests computing an approximation k to log10(&d) by
220 *
221 * k = (i - Bias)*0.301029995663981
222 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
223 *
224 * We want k to be too large rather than too small.
225 * The error in the first-order Taylor series approximation
226 * is in our favor, so we just round up the constant enough
227 * to compensate for any error in the multiplication of
228 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
229 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
230 * adding 1e-13 to the constant term more than suffices.
231 * Hence we adjust the constant term to 0.1760912590558.
232 * (We could get a more accurate k by invoking log10,
233 * but this is probably not worthwhile.)
234 */
235
236 i -= Bias;
237 #ifdef IBM
238 i <<= 2;
239 i += j;
240 #endif
241 #ifndef Sudden_Underflow
242 denorm = 0;
243 }
244 else {
245 /* d is denormalized */
246
247 i = bbits + be + (Bias + (P-1) - 1);
248 x = i > 32 ? word0(&d) << (64 - i) | word1(&d) >> (i - 32)
249 : word1(&d) << (32 - i);
250 dval(&d2) = x;
251 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
252 i -= (Bias + (P-1) - 1) + 1;
253 denorm = 1;
254 }
255 #endif
256 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
257 k = (int)ds;
258 if (ds < 0. && ds != k)
259 k--; /* want k = floor(ds) */
260 k_check = 1;
261 if (k >= 0 && k <= Ten_pmax) {
262 if (dval(&d) < tens[k])
263 k--;
264 k_check = 0;
265 }
266 j = bbits - i - 1;
267 if (j >= 0) {
268 b2 = 0;
269 s2 = j;
270 }
271 else {
272 b2 = -j;
273 s2 = 0;
274 }
275 if (k >= 0) {
276 b5 = 0;
277 s5 = k;
278 s2 += k;
279 }
280 else {
281 b2 -= k;
282 b5 = -k;
283 s5 = 0;
284 }
285 if (mode < 0 || mode > 9)
286 mode = 0;
287
288 #ifndef SET_INEXACT
289 #ifdef Check_FLT_ROUNDS
290 try_quick = Rounding == 1;
291 #else
292 try_quick = 1;
293 #endif
294 #endif /*SET_INEXACT*/
295
296 if (mode > 5) {
297 mode -= 4;
298 try_quick = 0;
299 }
300 leftright = 1;
301 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
302 /* silence erroneous "gcc -Wall" warning. */
303 switch(mode) {
304 case 0:
305 case 1:
306 i = 18;
307 ndigits = 0;
308 break;
309 case 2:
310 leftright = 0;
311 /* no break */
312 case 4:
313 if (ndigits <= 0)
314 ndigits = 1;
315 ilim = ilim1 = i = ndigits;
316 break;
317 case 3:
318 leftright = 0;
319 /* no break */
320 case 5:
321 i = ndigits + k + 1;
322 ilim = i;
323 ilim1 = i - 1;
324 if (i <= 0)
325 i = 1;
326 }
327 s = s0 = rv_alloc(i);
328
329 #ifdef Honor_FLT_ROUNDS
330 if (mode > 1 && Rounding != 1)
331 leftright = 0;
332 #endif
333
334 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
335
336 /* Try to get by with floating-point arithmetic. */
337
338 i = 0;
339 dval(&d2) = dval(&d);
340 k0 = k;
341 ilim0 = ilim;
342 ieps = 2; /* conservative */
343 if (k > 0) {
344 ds = tens[k&0xf];
345 j = k >> 4;
346 if (j & Bletch) {
347 /* prevent overflows */
348 j &= Bletch - 1;
349 dval(&d) /= bigtens[n_bigtens-1];
350 ieps++;
351 }
352 for(; j; j >>= 1, i++)
353 if (j & 1) {
354 ieps++;
355 ds *= bigtens[i];
356 }
357 dval(&d) /= ds;
358 }
359 else if (( j1 = -k )!=0) {
360 dval(&d) *= tens[j1 & 0xf];
361 for(j = j1 >> 4; j; j >>= 1, i++)
362 if (j & 1) {
363 ieps++;
364 dval(&d) *= bigtens[i];
365 }
366 }
367 if (k_check && dval(&d) < 1. && ilim > 0) {
368 if (ilim1 <= 0)
369 goto fast_failed;
370 ilim = ilim1;
371 k--;
372 dval(&d) *= 10.;
373 ieps++;
374 }
375 dval(&eps) = ieps*dval(&d) + 7.;
376 word0(&eps) -= (P-1)*Exp_msk1;
377 if (ilim == 0) {
378 S = mhi = 0;
379 dval(&d) -= 5.;
380 if (dval(&d) > dval(&eps))
381 goto one_digit;
382 if (dval(&d) < -dval(&eps))
383 goto no_digits;
384 goto fast_failed;
385 }
386 #ifndef No_leftright
387 if (leftright) {
388 /* Use Steele & White method of only
389 * generating digits needed.
390 */
391 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
392 for(i = 0;;) {
393 L = dval(&d);
394 dval(&d) -= L;
395 *s++ = '0' + (int)L;
396 if (dval(&d) < dval(&eps))
397 goto ret1;
398 if (1. - dval(&d) < dval(&eps))
399 goto bump_up;
400 if (++i >= ilim)
401 break;
402 dval(&eps) *= 10.;
403 dval(&d) *= 10.;
404 }
405 }
406 else {
407 #endif
408 /* Generate ilim digits, then fix them up. */
409 dval(&eps) *= tens[ilim-1];
410 for(i = 1;; i++, dval(&d) *= 10.) {
411 L = (Long)(dval(&d));
412 if (!(dval(&d) -= L))
413 ilim = i;
414 *s++ = '0' + (int)L;
415 if (i == ilim) {
416 if (dval(&d) > 0.5 + dval(&eps))
417 goto bump_up;
418 else if (dval(&d) < 0.5 - dval(&eps)) {
419 while(*--s == '0');
420 s++;
421 goto ret1;
422 }
423 break;
424 }
425 }
426 #ifndef No_leftright
427 }
428 #endif
429 fast_failed:
430 s = s0;
431 dval(&d) = dval(&d2);
432 k = k0;
433 ilim = ilim0;
434 }
435
436 /* Do we have a "small" integer? */
437
438 if (be >= 0 && k <= Int_max) {
439 /* Yes. */
440 ds = tens[k];
441 if (ndigits < 0 && ilim <= 0) {
442 S = mhi = 0;
443 if (ilim < 0 || dval(&d) <= 5*ds)
444 goto no_digits;
445 goto one_digit;
446 }
447 for(i = 1;; i++, dval(&d) *= 10.) {
448 L = (Long)(dval(&d) / ds);
449 dval(&d) -= L*ds;
450 #ifdef Check_FLT_ROUNDS
451 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
452 if (dval(&d) < 0) {
453 L--;
454 dval(&d) += ds;
455 }
456 #endif
457 *s++ = '0' + (int)L;
458 if (!dval(&d)) {
459 #ifdef SET_INEXACT
460 inexact = 0;
461 #endif
462 break;
463 }
464 if (i == ilim) {
465 #ifdef Honor_FLT_ROUNDS
466 if (mode > 1)
467 switch(Rounding) {
468 case 0: goto ret1;
469 case 2: goto bump_up;
470 }
471 #endif
472 dval(&d) += dval(&d);
473 if (dval(&d) > ds || (dval(&d) == ds && L & 1)) {
474 bump_up:
475 while(*--s == '9')
476 if (s == s0) {
477 k++;
478 *s = '0';
479 break;
480 }
481 ++*s++;
482 }
483 break;
484 }
485 }
486 goto ret1;
487 }
488
489 m2 = b2;
490 m5 = b5;
491 mhi = mlo = 0;
492 if (leftright) {
493 i =
494 #ifndef Sudden_Underflow
495 denorm ? be + (Bias + (P-1) - 1 + 1) :
496 #endif
497 #ifdef IBM
498 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
499 #else
500 1 + P - bbits;
501 #endif
502 b2 += i;
503 s2 += i;
504 mhi = i2b(1);
505 }
506 if (m2 > 0 && s2 > 0) {
507 i = m2 < s2 ? m2 : s2;
508 b2 -= i;
509 m2 -= i;
510 s2 -= i;
511 }
512 if (b5 > 0) {
513 if (leftright) {
514 if (m5 > 0) {
515 mhi = pow5mult(mhi, m5);
516 b1 = mult(mhi, b);
517 Bfree(b);
518 b = b1;
519 }
520 if (( j = b5 - m5 )!=0)
521 b = pow5mult(b, j);
522 }
523 else
524 b = pow5mult(b, b5);
525 }
526 S = i2b(1);
527 if (s5 > 0)
528 S = pow5mult(S, s5);
529
530 /* Check for special case that d is a normalized power of 2. */
531
532 spec_case = 0;
533 if ((mode < 2 || leftright)
534 #ifdef Honor_FLT_ROUNDS
535 && Rounding == 1
536 #endif
537 ) {
538 if (!word1(&d) && !(word0(&d) & Bndry_mask)
539 #ifndef Sudden_Underflow
540 && word0(&d) & (Exp_mask & ~Exp_msk1)
541 #endif
542 ) {
543 /* The special case */
544 b2 += Log2P;
545 s2 += Log2P;
546 spec_case = 1;
547 }
548 }
549
550 /* Arrange for convenient computation of quotients:
551 * shift left if necessary so divisor has 4 leading 0 bits.
552 *
553 * Perhaps we should just compute leading 28 bits of S once
554 * and for all and pass them and a shift to quorem, so it
555 * can do shifts and ors to compute the numerator for q.
556 */
557 #ifdef Pack_32
558 if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0)
559 i = 32 - i;
560 #else
561 if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0)
562 i = 16 - i;
563 #endif
564 if (i > 4) {
565 i -= 4;
566 b2 += i;
567 m2 += i;
568 s2 += i;
569 }
570 else if (i < 4) {
571 i += 28;
572 b2 += i;
573 m2 += i;
574 s2 += i;
575 }
576 if (b2 > 0)
577 b = lshift(b, b2);
578 if (s2 > 0)
579 S = lshift(S, s2);
580 if (k_check) {
581 if (cmp(b,S) < 0) {
582 k--;
583 b = multadd(b, 10, 0); /* we botched the k estimate */
584 if (leftright)
585 mhi = multadd(mhi, 10, 0);
586 ilim = ilim1;
587 }
588 }
589 if (ilim <= 0 && (mode == 3 || mode == 5)) {
590 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
591 /* no digits, fcvt style */
592 no_digits:
593 k = -1 - ndigits;
594 goto ret;
595 }
596 one_digit:
597 *s++ = '1';
598 k++;
599 goto ret;
600 }
601 if (leftright) {
602 if (m2 > 0)
603 mhi = lshift(mhi, m2);
604
605 /* Compute mlo -- check for special case
606 * that d is a normalized power of 2.
607 */
608
609 mlo = mhi;
610 if (spec_case) {
611 mhi = Balloc(mhi->k);
612 Bcopy(mhi, mlo);
613 mhi = lshift(mhi, Log2P);
614 }
615
616 for(i = 1;;i++) {
617 dig = quorem(b,S) + '0';
618 /* Do we yet have the shortest decimal string
619 * that will round to d?
620 */
621 j = cmp(b, mlo);
622 delta = diff(S, mhi);
623 j1 = delta->sign ? 1 : cmp(b, delta);
624 Bfree(delta);
625 #ifndef ROUND_BIASED
626 if (j1 == 0 && mode != 1 && !(word1(&d) & 1)
627 #ifdef Honor_FLT_ROUNDS
628 && Rounding >= 1
629 #endif
630 ) {
631 if (dig == '9')
632 goto round_9_up;
633 if (j > 0)
634 dig++;
635 #ifdef SET_INEXACT
636 else if (!b->x[0] && b->wds <= 1)
637 inexact = 0;
638 #endif
639 *s++ = dig;
640 goto ret;
641 }
642 #endif
643 if (j < 0 || (j == 0 && mode != 1
644 #ifndef ROUND_BIASED
645 && !(word1(&d) & 1)
646 #endif
647 )) {
648 if (!b->x[0] && b->wds <= 1) {
649 #ifdef SET_INEXACT
650 inexact = 0;
651 #endif
652 goto accept_dig;
653 }
654 #ifdef Honor_FLT_ROUNDS
655 if (mode > 1)
656 switch(Rounding) {
657 case 0: goto accept_dig;
658 case 2: goto keep_dig;
659 }
660 #endif /*Honor_FLT_ROUNDS*/
661 if (j1 > 0) {
662 b = lshift(b, 1);
663 j1 = cmp(b, S);
664 if ((j1 > 0 || (j1 == 0 && dig & 1))
665 && dig++ == '9')
666 goto round_9_up;
667 }
668 accept_dig:
669 *s++ = dig;
670 goto ret;
671 }
672 if (j1 > 0) {
673 #ifdef Honor_FLT_ROUNDS
674 if (!Rounding)
675 goto accept_dig;
676 #endif
677 if (dig == '9') { /* possible if i == 1 */
678 round_9_up:
679 *s++ = '9';
680 goto roundoff;
681 }
682 *s++ = dig + 1;
683 goto ret;
684 }
685 #ifdef Honor_FLT_ROUNDS
686 keep_dig:
687 #endif
688 *s++ = dig;
689 if (i == ilim)
690 break;
691 b = multadd(b, 10, 0);
692 if (mlo == mhi)
693 mlo = mhi = multadd(mhi, 10, 0);
694 else {
695 mlo = multadd(mlo, 10, 0);
696 mhi = multadd(mhi, 10, 0);
697 }
698 }
699 }
700 else
701 for(i = 1;; i++) {
702 *s++ = dig = quorem(b,S) + '0';
703 if (!b->x[0] && b->wds <= 1) {
704 #ifdef SET_INEXACT
705 inexact = 0;
706 #endif
707 goto ret;
708 }
709 if (i >= ilim)
710 break;
711 b = multadd(b, 10, 0);
712 }
713
714 /* Round off last digit */
715
716 #ifdef Honor_FLT_ROUNDS
717 switch(Rounding) {
718 case 0: goto trimzeros;
719 case 2: goto roundoff;
720 }
721 #endif
722 b = lshift(b, 1);
723 j = cmp(b, S);
724 if (j > 0 || (j == 0 && dig & 1)) {
725 roundoff:
726 while(*--s == '9')
727 if (s == s0) {
728 k++;
729 *s++ = '1';
730 goto ret;
731 }
732 ++*s++;
733 }
734 else {
735 #ifdef Honor_FLT_ROUNDS
736 trimzeros:
737 #endif
738 while(*--s == '0');
739 s++;
740 }
741 ret:
742 Bfree(S);
743 if (mhi) {
744 if (mlo && mlo != mhi)
745 Bfree(mlo);
746 Bfree(mhi);
747 }
748 ret1:
749 #ifdef SET_INEXACT
750 if (inexact) {
751 if (!oldinexact) {
752 word0(&d) = Exp_1 + (70 << Exp_shift);
753 word1(&d) = 0;
754 dval(&d) += 1.;
755 }
756 }
757 else if (!oldinexact)
758 clear_inexact();
759 #endif
760 Bfree(b);
761 *s = 0;
762 *decpt = k + 1;
763 if (rve)
764 *rve = s;
765 return s0;
766 }