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