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