8322 nl: misleading-indentation
[unleashed/tickless.git] / usr / src / lib / libbc / libc / gen / common / double_decim.c
blob53be17331b49d6d464487024af92aaf3ba05527d
1 /*
2 * CDDL HEADER START
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License, Version 1.0 only
6 * (the "License"). You may not use this file except in compliance
7 * with the License.
9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
10 * or http://www.opensolaris.org/os/licensing.
11 * See the License for the specific language governing permissions
12 * and limitations under the License.
14 * When distributing Covered Code, include this CDDL HEADER in each
15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
16 * If applicable, add the following below this CDDL HEADER, with the
17 * fields enclosed by brackets "[]" replaced with your own identifying
18 * information: Portions Copyright [yyyy] [name of copyright owner]
20 * CDDL HEADER END
23 * Copyright 1988 Sun Microsystems, Inc. All rights reserved.
24 * Use is subject to license terms.
27 #pragma ident "%Z%%M% %I% %E% SMI"
29 /* Conversion between binary and decimal floating point. */
31 #include "base_conversion.h"
33 /* PRIVATE FUNCTIONS */
36 * Rounds decimal record *pd according to modes in *pm, recording exceptions
37 * for inexact or overflow in *ps. round is the round digit and sticky is 0
38 * or non-zero to indicate exact or inexact. pd->ndigits is expected to be
39 * correctly set.
41 void
42 decimal_round(decimal_mode *pm, decimal_record *pd, fp_exception_field_type *ps,
43 char round, unsigned sticky)
45 int lsd, i;
47 if ((round == '0') && (sticky == 0)) { /* Exact. */
48 goto done;
50 *ps |= 1 << fp_inexact;
52 switch (pm->rd) {
53 case fp_nearest:
54 if (round < '5')
55 goto done;
56 if (round > '5')
57 goto roundup;
58 if (sticky != 0)
59 goto roundup;
60 /* Now in ambiguous case; round up if lsd is odd. */
61 if (pd->ndigits <= 0)
62 goto done; /* Presumed 0. */
63 lsd = pd->ds[pd->ndigits - 1] - '0';
64 if ((lsd % 2) == 0)
65 goto done;
66 goto roundup;
67 case fp_positive:
68 if (pd->sign != 0)
69 goto done;
70 goto roundup;
71 case fp_negative:
72 if (pd->sign == 0)
73 goto done;
74 goto roundup;
75 case fp_tozero:
76 goto done;
78 roundup:
79 for (i = (pd->ndigits - 1); (pd->ds[i] == '9') && (i >= 0); i--)
80 pd->ds[i] = '0';
81 if (i >= 0)
82 (pd->ds[i])++;
83 else { /* Rounding carry out has occurred. */
84 pd->ds[0] = '1';
85 if (pm->df == floating_form) { /* For E format, simply
86 * adjust exponent. */
87 pd->exponent++;
88 } else { /* For F format, increase length of string. */
89 if (pd->ndigits > 0)
90 pd->ds[pd->ndigits] = '0';
91 pd->ndigits++;
94 goto ret;
95 done:
96 if (pd->ndigits <= 0) { /* Create zero string. */
97 pd->ds[0] = '0';
98 pd->ndigits = 1;
100 ret:
101 pd->ds[pd->ndigits] = 0;/* Terminate string. */
102 return;
106 * Converts an unpacked integer value *pu into a decimal string in *ds, of
107 * length returned in *ndigs. Inexactness is indicated by setting
108 * ds[ndigs-1] odd.
110 * Arguments
111 * pu: Input unpacked integer value input.
112 * nsig: Input number of significant digits required.
113 * ds: Output decimal integer string output
114 * must be large enough.
115 * nzeros: Output number of implicit trailing zeros
116 * produced.
117 * ndigs: Output number of explicit digits produced
118 * in ds.
120 void
121 binary_to_decimal_integer(unpacked *pu, unsigned nsig, char ds[],
122 unsigned *nzeros, unsigned *ndigs)
125 _big_float *pd, b, d;
126 int e, i, is;
127 _BIG_FLOAT_DIGIT stickyshift;
128 char s[4];
130 b.bsize = _BIG_FLOAT_SIZE; /* Initialize sizes of big floats. */
131 d.bsize = _BIG_FLOAT_SIZE;
132 _unpacked_to_big_float(pu, &b, &e);
133 if (e < 0) {
134 _right_shift_base_two(&b, (short unsigned) -e, &stickyshift);
135 #ifdef DEBUG
136 assert(stickyshift == 0);
137 #endif
139 _big_binary_to_big_decimal(&b, &d);
140 if (e <= 0)
141 pd = &d;
142 else {
143 _big_float_times_power(&d, 2, e, (int) nsig, &pd);
144 switch ((unsigned int)pd) {
145 case ((unsigned int)BIG_FLOAT_TIMES_TOOBIG):
147 char bcastring[80];
149 (void) sprintf(bcastring, " binary exponent %d ", e);
150 _base_conversion_abort(ERANGE, bcastring);
151 break;
153 case ((unsigned int)BIG_FLOAT_TIMES_NOMEM):
155 char bcastring[80];
157 (void) sprintf(bcastring, " binary exponent %d ", e);
158 _base_conversion_abort(ENOMEM, bcastring);
159 break;
161 default:
162 #ifdef DEBUG
163 if (pd != &d)
164 (void) printf(" large binary exponent %d needs heap buffer \n", e);
165 printf(" product ");
166 _display_big_float(pd, 10);
167 #endif
168 break;
171 _fourdigitsquick((short unsigned) pd->bsignificand[pd->blength - 1], s);
172 for (i = 0; s[i] == '0'; i++); /* Find first non-zero digit. */
173 for (is = 0; i <= 3;)
174 ds[is++] = s[i++]; /* Copy subsequent digits. */
176 for (i = (pd->blength - 2); i >= 0; i--) { /* Convert powers of
177 * 10**4 to decimal
178 * digits. */
179 _fourdigitsquick((short unsigned) pd->bsignificand[i], &(ds[is]));
180 is += 4;
183 ds[is] = 0;
184 *ndigs = is;
185 *nzeros = pd->bexponent;
186 if (pd != &d)
187 _free_big_float(pd);
189 #ifdef DEBUG
190 printf(" binary to decimal integer result %s * 10**%d \n", ds, *nzeros);
191 #endif
195 * Converts an unpacked fraction value *pu into a decimal string consisting
196 * of a) an implicit '.' b) *nzeros implicit leading zeros c) *ndigs explicit
197 * digits in ds ds contains at least nsig significant digits. nzeros + *
198 * *ndigs is at least nfrac digits after the point. Inexactness is indicated
199 * by sticking to the lsb.
201 * Arguments
203 * pu: Input unpacked fraction value output < 1
204 * in magnitude.
205 * nsig: Input number of significant digits
206 * required.
207 * nfrac: Input number of digits after point
208 * required.
209 * ds: Output decimal integer string output -
210 * must be large enough.
211 * nzeros: Output number of implicit leading zeros
212 * produced.
213 * ndigs: Output number of explicit digits produced
214 * in ds.
216 void
217 binary_to_decimal_fraction(unpacked *pu, unsigned nsig, unsigned nfrac,
218 char ds[], int *nzeros, int *ndigs)
220 _big_float *pb, b, d;
221 int e, i, j, is, excess;
222 char s[4];
223 int tensig, tenpower;
224 _BIG_FLOAT_DIGIT stickyshift;
226 *nzeros = 0;
227 if (pu->fpclass == fp_zero) { /* Exact zero. */
228 for (i = 0; i <= nfrac; i++)
229 ds[i] = '0';
230 for (; i <= nsig; i++)
231 ds[i] = '0';
232 *ndigs = i;
233 return;
235 b.bsize = _BIG_FLOAT_SIZE; /* Initialize sizes of big floats. */
236 d.bsize = _BIG_FLOAT_SIZE;
237 _unpacked_to_big_float(pu, &b, &e);
239 * e < 0 always
241 b.bexponent = e;
242 tenpower = nsig + (int) (((17 - e - 16 * b.blength) * (unsigned long) 19729) >> 16);
243 if (tenpower < nfrac)
244 tenpower = nfrac;
245 tensig = nfrac;
246 if (nsig > tensig)
247 tensig = nsig;
248 tensig = 1 + (((tensig + 2) * 217706) >> 16);
249 tensig = -tensig;
251 #ifdef DEBUG
252 printf(" binary to decimal fraction exponent 2**%d \n", e);
253 printf(" binary to decimal fraction nsig %d nfrac %d tenpower %d tensig %d \n", nsig, nfrac, tenpower, tensig);
254 #endif
255 _big_float_times_power(&b, 10, tenpower, tensig, &pb);
256 switch ((unsigned int)pb) {
257 case ((unsigned int)BIG_FLOAT_TIMES_TOOBIG):
259 char bcastring[80];
261 (void) sprintf(bcastring, " decimal exponent %d ", tenpower);
262 _base_conversion_abort(ERANGE, bcastring);
263 break;
265 case ((unsigned int)BIG_FLOAT_TIMES_NOMEM):
267 char bcastring[80];
269 (void) sprintf(bcastring, " decimal exponent %d ", tenpower);
270 _base_conversion_abort(ENOMEM, bcastring);
271 break;
273 default:
274 #ifdef DEBUG
275 if (pb != &b)
276 printf(" large decimal exponent %d needs heap buffer \n", tenpower);
277 printf(" product ");
278 _display_big_float(pb, 2);
279 #endif
280 break;
283 if (pb->bexponent <= -16) {
284 /* Have computed appropriate decimal part; now toss fraction. */
285 excess = (-pb->bexponent) / 16;
286 #ifdef DEBUG
287 printf(" discard %d excess fraction bits \n", 16 * excess);
288 #endif
289 for (i = 0; (i < excess) && (pb->bsignificand[i] == 0); i++);
290 if (i < excess)
291 pb->bsignificand[excess] |= 1; /* Sticky bit for
292 * discarded fraction. */
293 for (i = excess; i < pb->blength; i++)
294 pb->bsignificand[i - excess] = pb->bsignificand[i];
296 pb->blength -= excess;
297 pb->bexponent += 16 * excess;
299 if (pb->bexponent < 0) {
300 _right_shift_base_two(pb, (short unsigned) -pb->bexponent, &stickyshift);
301 if (stickyshift != 0)
302 pb->bsignificand[0] |= 1; /* Stick to lsb. */
304 _big_binary_to_big_decimal(pb, &d);
306 i = d.blength - 1;
307 while (d.bsignificand[i] == 0)
308 i--;
309 _fourdigitsquick((short unsigned) d.bsignificand[i], s);
310 for (j = 0; s[j] == '0'; j++); /* Find first non-zero digit. */
311 for (is = 0; j <= 3;)
312 ds[is++] = s[j++]; /* Copy subsequent digits. */
314 for (i--; i >= 0; i--) {/* Convert powers of 10**4 to decimal digits. */
315 _fourdigitsquick((short unsigned) d.bsignificand[i], &(ds[is]));
316 is += 4;
319 ds[is] = 0;
320 *ndigs = is;
321 #ifdef DEBUG
322 assert(tenpower >= is);
323 #endif
324 *nzeros = tenpower - is;/* There were supposed to be tenpower leading
325 * digits, and is were found. */
327 if (pb != &b)
328 _free_big_float(pb);
330 #ifdef DEBUG
331 printf(" binary to decimal fraction result .%s * 10**%d \n", ds, -(*nzeros));
332 #endif
336 void
337 _unpacked_to_decimal(unpacked *px, decimal_mode *pm, decimal_record *pd,
338 fp_exception_field_type *ps)
340 unpacked fx, ix;
341 unsigned fmask, imask;
342 int i, intdigs, fracdigs, fraczeros, fracsigs, ids, idsbound, lzbound;
343 unsigned nsig, nfrac, intzeros, intsigs;
344 char is[_INTEGER_SIZE], fs[DECIMAL_STRING_LENGTH];
345 char round = '0';
346 unsigned sticky = 0;
348 pd->sign = px->sign;
349 pd->fpclass = px->fpclass;
350 if ((px->fpclass != fp_normal) && (px->fpclass != fp_subnormal))
351 return;
352 if ((pm->ndigits >= DECIMAL_STRING_LENGTH) ||
353 ((pm->df == floating_form) && (pm->ndigits < 1))) { /* Gross overflow or bad
354 * spec. */
355 overflow:
356 *ps |= 1 << fp_overflow;
357 return;
359 /* Divide x up into integer part ix and fraction part fx. */
361 ix = *px;
362 fx = ix;
363 if (ix.exponent <= -1) {/* All fraction. */
364 ix.fpclass = fp_zero;
365 } else if (ix.exponent >= 159) { /* All integer. */
366 fx.fpclass = fp_zero;
367 } else if ((ix.exponent % 32) == 31) { /* Integer/fraction boundary
368 * is conveniently on a word
369 * boundary. */
370 imask = (ix.exponent + 1) / 32; /* Words 0..imask-1 are
371 * integer; imask..SIZE are
372 * fraction. */
373 for (i = 0; i < imask; i++)
374 fx.significand[i] = 0;
375 for (; i < UNPACKED_SIZE; i++)
376 ix.significand[i] = 0;
377 _fp_normalize(&fx);
378 } else { /* Integer/fraction boundary falls in the
379 * middle of a word. */
380 imask = (ix.exponent + 1) / 32; /* Words 0..imask-1 are
381 * integer; imask is integer
382 * and fraction ;
383 * imask+1..SIZE are
384 * fraction. */
385 for (i = 0; i < imask; i++)
386 fx.significand[i] = 0;
387 fmask = (1 << (31 - (ix.exponent % 32))) - 1;
388 fx.significand[imask] &= fmask;
389 ix.significand[imask] &= ~fmask;
390 for (i = (imask + 1); i < UNPACKED_SIZE; i++)
391 ix.significand[i] = 0;
392 _fp_normalize(&fx);
394 if (ix.fpclass != fp_zero) { /* Compute integer part of result. */
395 if (pm->df == floating_form)
396 nsig = pm->ndigits + 1; /* Significant digits wanted
397 * for E format, plus one for
398 * rounding. */
399 else
400 nsig = _INTEGER_SIZE; /* Significant digits wanted
401 * for F format == all. */
403 binary_to_decimal_integer(&ix, nsig, is, &intzeros, &intsigs);
404 } else {
405 intsigs = 0;
406 intzeros = 0;
408 intdigs = intsigs + intzeros;
409 fracdigs = 0;
410 if (((pm->df == fixed_form) && (pm->ndigits >= 0)) ||
411 ((pm->df == floating_form) && ((pm->ndigits + 1) > intdigs))) { /* Need to compute
412 * fraction part. */
413 if (pm->df == floating_form) { /* Need more significant
414 * digits. */
415 nsig = pm->ndigits + 2 - intdigs; /* Add two for rounding,
416 * sticky. */
417 if (nsig > DECIMAL_STRING_LENGTH)
418 nsig = DECIMAL_STRING_LENGTH;
419 nfrac = 1;
420 } else { /* Need fraction digits. */
421 nsig = 0;
422 nfrac = pm->ndigits + 2; /* Add two for rounding,
423 * sticky. */
424 if (nfrac > DECIMAL_STRING_LENGTH)
425 nfrac = DECIMAL_STRING_LENGTH;
427 binary_to_decimal_fraction(&fx, nsig, nfrac, fs, &fraczeros, &fracsigs);
428 fracdigs = fraczeros + fracsigs;
430 if (pm->df == floating_form) { /* Combine integer and fraction for E
431 * format. */
432 idsbound = intsigs;
433 if (idsbound > pm->ndigits)
434 idsbound = pm->ndigits;
435 for (ids = 0; ids < idsbound; ids++)
436 pd->ds[ids] = is[ids];
437 /* Put integer into output string. */
438 idsbound = intsigs + intzeros;
439 if (idsbound > pm->ndigits)
440 idsbound = pm->ndigits;
441 for (; ids < idsbound; ids++)
442 pd->ds[ids] = '0';
443 if (ids == pm->ndigits) { /* Integer part had enough
444 * significant digits. */
445 pd->ndigits = ids;
446 pd->exponent = intdigs - ids;
447 if (ids < intdigs) { /* Gather rounding info. */
448 if (ids < intsigs)
449 round = is[ids++];
450 else
451 round = '0';
452 for (; (is[ids] == '0') && (ids < intsigs); ids++);
453 if (ids < intsigs)
454 sticky = 1;
455 if (fx.fpclass != fp_zero)
456 sticky = 1;
457 } else {/* Integer part is exact - round from
458 * fraction. */
459 if (fx.fpclass != fp_zero) {
460 int stickystart;
461 /* Fraction non-zero. */
462 if (fraczeros > 0) { /* Round digit is zero. */
463 round = '0';
464 stickystart = 0; /* Stickies start with
465 * fs[0]. */
466 } else { /* Round digit is fs[0]. */
467 round = fs[0];
468 stickystart = 1; /* Stickies start with
469 * fs[1]. */
471 if (sticky == 0) { /* Search for sticky
472 * bits. */
473 for (ids = stickystart; (fs[ids] == '0') && (ids < fracdigs); ids++);
474 if (ids < fracdigs)
475 sticky = 1;
479 } else { /* Need more significant digits from fraction
480 * part. */
481 idsbound = pm->ndigits - ids;
482 if (ids == 0) { /* No integer part - find first
483 * significant digit. */
484 for (i = 0; fs[i] == '0'; i++);
485 idsbound = i + idsbound + fraczeros;
486 i += fraczeros; /* Point i at first
487 * significant digit. */
488 } else
489 i = 0;
490 if (idsbound > fracdigs)
491 idsbound = fracdigs;
492 pd->exponent = -idsbound;
494 if (fraczeros < idsbound) /* Compute number of
495 * leading zeros
496 * required. */
497 lzbound = fraczeros;
498 else
499 lzbound = idsbound;
500 for (; (i < lzbound); i++)
501 pd->ds[ids++] = '0';
502 for (; (i < idsbound); i++)
503 pd->ds[ids++] = fs[i - fraczeros];
504 i -= fraczeros; /* Don't worry about leading zeros
505 * from now on, we're just rounding */
506 if (i < fracsigs) { /* Gather rounding info. */
507 if (i < 0)
508 round = '0';
509 else
510 round = fs[i];
511 i++;
512 if (sticky == 0) { /* Find out if remainder
513 * is exact. */
514 if (i < 0)
515 i = 0;
516 for (; (fs[i] == '0') && (i < fracsigs); i++);
517 if (i < fracsigs)
518 sticky = 1;
520 } else {/* Fraction part is exact - add zero digits
521 * if required. */
522 for (; ids < pm->ndigits; ids++)
523 pd->ds[ids] = '0';
525 pd->ndigits = ids;
527 decimal_round(pm, pd, ps, round, sticky);
528 } else { /* Combine integer and fraction for F format. */
529 if (pm->ndigits >= 0) { /* Normal F format. */
530 if ((intdigs + pm->ndigits) >= DECIMAL_STRING_LENGTH)
531 goto overflow;
532 for (ids = 0; ids < intsigs; ids++)
533 pd->ds[ids] = is[ids];
534 for (; ids < intdigs; ids++)
535 pd->ds[ids] = '0';
536 /* Copy integer digits. */
537 idsbound = fracdigs;
538 if (idsbound > pm->ndigits)
539 idsbound = pm->ndigits;
540 if (fraczeros < idsbound) /* Compute number of
541 * leading zeros
542 * required. */
543 lzbound = fraczeros;
544 else
545 lzbound = idsbound;
546 for (i = 0; (i < lzbound); i++)
547 pd->ds[ids++] = '0';
548 for (; (i < idsbound); i++)
549 pd->ds[ids++] = fs[i - fraczeros]; /* Copy fraction digits. */
550 for (; i < pm->ndigits; i++)
551 pd->ds[ids++] = '0';
552 /* Copy trailing zeros if necessary. */
553 pd->ndigits = ids;
554 pd->exponent = intdigs - ids;
555 i -= fraczeros; /* Don't worry about leading zeros
556 * from now on, we're just rounding */
557 if (i < fracsigs) { /* Gather rounding info. */
558 if (i < 0)
559 round = '0';
560 else
561 round = fs[i];
562 i++;
563 if (sticky == 0) { /* Find out if remainder
564 * is exact. */
565 if (i < 0)
566 i = 0;
567 for (; (fs[i] == '0') && (i < fracsigs); i++);
568 if (i < fracsigs)
569 sticky = 1;
572 decimal_round(pm, pd, ps, round, sticky);
573 } else { /* Bizarre F format - round to left of point. */
574 int roundpos = -pm->ndigits;
576 if (intdigs >= DECIMAL_STRING_LENGTH)
577 goto overflow;
578 if (roundpos >= DECIMAL_STRING_LENGTH)
579 goto overflow;
580 if (intdigs <= roundpos) { /* Not enough integer
581 * digits. */
582 if (intdigs == roundpos) {
583 round = is[0];
584 i = 1;
585 } else {
586 round = '0';
587 i = 0;
589 for (; (is[i] == '0') && (i < intsigs); i++);
590 /* Search for sticky bits. */
591 if (i < intsigs)
592 sticky = 1;
593 pd->ndigits = 0;
594 } else {/* Some integer digits do not get rounded
595 * away. */
596 #ifdef _NO_GOOD
597 for (ids = 0; ids < (intsigs - roundpos); ids++)
598 pd->ds[ids] = is[ids];
599 for (ids = 0; ids < (intdigs - roundpos); ids++)
600 pd->ds[ids] = '0';
601 #else
603 int ncopy = intsigs - roundpos;
604 if (ncopy > 0) {
605 /* Copy integer digits. */
606 (void) memcpy(&(pd->ds[0]), &(is[0]), ncopy);
607 ids = ncopy;
611 int ncopy = intdigs - roundpos - ids ;
612 if (ncopy > 0) {
613 (void) memset(&(pd->ds[ids]), '0', ncopy);
614 ids += ncopy;
617 #endif /* _NO_GOOD */
618 /* Copy integer digits. */
619 pd->ndigits = ids;
620 if (ids < intsigs) { /* Inexact. */
621 round = is[ids++];
622 for (; (is[ids] == '0') && (ids < intsigs); ids++);
623 /* Search for non-zero digits. */
624 if (ids < intsigs)
625 sticky = 1;
628 if (fx.fpclass != fp_zero)
629 sticky = 1;
630 decimal_round(pm, pd, ps, round, sticky);
631 for (i = pd->ndigits; i < (pd->ndigits + roundpos); i++)
632 pd->ds[i] = '0'; /* Blank out rounded
633 * away digits. */
634 pd->exponent = 0;
635 pd->ndigits = i;
636 pd->ds[i] = 0; /* Terminate string. */
641 void
642 double_to_decimal(double *px, decimal_mode *pm, decimal_record *pd,
643 fp_exception_field_type *ps)
645 double_equivalence kluge;
646 unpacked u;
648 *ps = 0; /* Initialize *ps. */
649 kluge.x = *px;
650 pd->sign = kluge.f.msw.sign;
651 pd->fpclass = _class_double(px);
652 switch (pd->fpclass) {
653 case fp_zero:
654 break;
655 case fp_infinity:
656 break;
657 case fp_quiet:
658 break;
659 case fp_signaling:
660 break;
661 default:
662 _unpack_double(&u, &kluge.x);
663 _unpacked_to_decimal(&u, pm, pd, ps);
667 void
668 quadruple_to_decimal(quadruple *px, decimal_mode *pm, decimal_record *pd,
669 fp_exception_field_type *ps)
671 quadruple_equivalence kluge;
672 unpacked u;
673 int i;
675 *ps = 0; /* Initialize *ps - no exceptions. */
676 for (i = 0; i < 4; i++)
677 #ifdef __STDC__
678 kluge.x = *px;
679 #else
680 kluge.x.u[i] = px->u[i];
681 #endif
682 pd->sign = kluge.f.msw.sign;
683 pd->fpclass = _class_quadruple(px);
684 switch (pd->fpclass) {
685 case fp_zero:
686 break;
687 case fp_infinity:
688 break;
689 case fp_quiet:
690 break;
691 case fp_signaling:
692 break;
693 default:
694 _unpack_quadruple(&u, px);
695 _unpacked_to_decimal(&u, pm, pd, ps);