xref: /illumos-gate/usr/src/lib/libc/port/fp/decimal_bin.c (revision 7257d1b4)
1 /*
2  * CDDL HEADER START
3  *
4  * The contents of this file are subject to the terms of the
5  * Common Development and Distribution License (the "License").
6  * You may not use this file except in compliance with the License.
7  *
8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9  * or http://www.opensolaris.org/os/licensing.
10  * See the License for the specific language governing permissions
11  * and limitations under the License.
12  *
13  * When distributing Covered Code, include this CDDL HEADER in each
14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15  * If applicable, add the following below this CDDL HEADER, with the
16  * fields enclosed by brackets "[]" replaced with your own identifying
17  * information: Portions Copyright [yyyy] [name of copyright owner]
18  *
19  * CDDL HEADER END
20  */
21 
22 /*
23  * Copyright 2008 Sun Microsystems, Inc.  All rights reserved.
24  * Use is subject to license terms.
25  */
26 
27 #pragma ident	"%Z%%M%	%I%	%E% SMI"
28 
29 /*
30  * Conversion from decimal to binary floating point
31  */
32 
33 #include "lint.h"
34 #include <stdlib.h>
35 #include "base_conversion.h"
36 
37 /*
38  * Convert the integer part of a nonzero base-10^4 _big_float *pd
39  * to base 2^16 in **ppb.  The converted value is accurate to nsig
40  * significant bits.  On exit, *sticky is nonzero if *pd had a
41  * nonzero fractional part.  If pd->exponent > 0 and **ppb is not
42  * large enough to hold the final converted value (i.e., the con-
43  * verted significand scaled by 10^pd->exponent), then on exit,
44  * *ppb will point to a newly allocated _big_float, which must be
45  * freed by the caller.  (The number of significant bits we need
46  * should fit in pb, but __big_float_times_power may allocate new
47  * storage anyway because the exact product could require more than
48  * 16000 bits.)
49  *
50  * This routine does not check that **ppb is large enough to hold
51  * the result of converting the significand of *pd.
52  */
53 static void
54 __big_decimal_to_big_binary(_big_float *pd, int nsig, _big_float **ppb,
55     int *sticky)
56 {
57 	_big_float	*pb;
58 	int		i, j, len, s;
59 	unsigned int	carry;
60 
61 	pb = *ppb;
62 
63 	/* convert pd a digit at a time, most significant first */
64 	if (pd->bexponent + ((pd->blength - 1) << 2) >= 0) {
65 		pb->bsignificand[0] = pd->bsignificand[pd->blength - 1];
66 		len = 1;
67 		for (i = pd->blength - 2; i >= 0 &&
68 		    pd->bexponent + (i << 2) >= 0; i--) {
69 			/* multiply pb by 10^4 and add next digit */
70 			carry = pd->bsignificand[i];
71 			for (j = 0; j < len; j++) {
72 				carry += (unsigned int)pb->bsignificand[j]
73 				    * 10000;
74 				pb->bsignificand[j] = carry & 0xffff;
75 				carry >>= 16;
76 			}
77 			if (carry)
78 				pb->bsignificand[j++] = carry;
79 			len = j;
80 		}
81 	} else {
82 		i = pd->blength - 1;
83 		len = 0;
84 	}
85 
86 	/* convert any partial digit */
87 	if (i >= 0 && pd->bexponent + (i << 2) > -4) {
88 		s = pd->bexponent + (i << 2) + 4;
89 		/* multiply pb by 10^s and add partial digit */
90 		carry = pd->bsignificand[i];
91 		if (s == 1) {
92 			s = carry % 1000;
93 			carry = carry / 1000;
94 			for (j = 0; j < len; j++) {
95 				carry += (unsigned int)pb->bsignificand[j]
96 				    * 10;
97 				pb->bsignificand[j] = carry & 0xffff;
98 				carry >>= 16;
99 			}
100 		} else if (s == 2) {
101 			s = carry % 100;
102 			carry = carry / 100;
103 			for (j = 0; j < len; j++) {
104 				carry += (unsigned int)pb->bsignificand[j]
105 				    * 100;
106 				pb->bsignificand[j] = carry & 0xffff;
107 				carry >>= 16;
108 			}
109 		} else {
110 			s = carry % 10;
111 			carry = carry / 10;
112 			for (j = 0; j < len; j++) {
113 				carry += (unsigned int)pb->bsignificand[j]
114 				    * 1000;
115 				pb->bsignificand[j] = carry & 0xffff;
116 				carry >>= 16;
117 			}
118 		}
119 		if (carry)
120 			pb->bsignificand[j++] = carry;
121 		len = j;
122 		i--;
123 	} else {
124 		s = 0;
125 	}
126 
127 	pb->blength = len;
128 	pb->bexponent = 0;
129 
130 	/* continue accumulating sticky flag */
131 	while (i >= 0)
132 		s |= pd->bsignificand[i--];
133 	*sticky = s;
134 
135 	if (pd->bexponent > 0) {
136 		/* scale pb by 10^pd->exponent */
137 		__big_float_times_power(pb, 10, pd->bexponent, nsig, ppb);
138 	}
139 }
140 
141 /*
142  * Convert the decimal_record *pd to an unpacked datum *px accurately
143  * enough that *px can be rounded correctly to sigbits significant bits.
144  * (We may assume sigbits <= 113.)
145  */
146 static void
147 __decimal_to_unpacked(unpacked *px, decimal_record *pd, int sigbits)
148 {
149 	_big_float	d, b, *pbd, *pbb;
150 	char		*ds;
151 	int		ids, i, ix, exp, ndigs;
152 	int		sticky, powtwo, sigdigits;
153 
154 	px->sign = pd->sign;
155 	px->fpclass = pd->fpclass;
156 	ds = pd->ds;
157 	ndigs = pd->ndigits;
158 	exp = pd->exponent;
159 
160 	/* remove trailing zeroes */
161 	while (ndigs > 0 && ds[ndigs - 1] == '0') {
162 		exp++;
163 		ndigs--;
164 	}
165 	if (ndigs < 1) {
166 		/* nothing left */
167 		px->fpclass = fp_zero;
168 		return;
169 	}
170 
171 	/* convert remaining digits to a base-10^4 _big_float */
172 	d.bsize = _BIG_FLOAT_SIZE;
173 	d.bexponent = exp;
174 	d.blength = (ndigs + 3) >> 2;
175 	i = d.blength - 1;
176 	ids = ndigs - (d.blength << 2);
177 	switch (ids) {
178 	case -1:
179 		d.bsignificand[i] = 100 * ds[ids + 1] +
180 		    10 * ds[ids + 2] + ds[ids + 3] - 111 * '0';
181 		i--;
182 		ids += 4;
183 		break;
184 
185 	case -2:
186 		d.bsignificand[i] = 10 * ds[ids + 2] + ds[ids + 3] - 11 * '0';
187 		i--;
188 		ids += 4;
189 		break;
190 
191 	case -3:
192 		d.bsignificand[i] = ds[ids + 3] - '0';
193 		i--;
194 		ids += 4;
195 		break;
196 	}
197 	while (i >= 0) {
198 		d.bsignificand[i] = 1000 * ds[ids] + 100 * ds[ids + 1] +
199 		    10 * ds[ids + 2] + ds[ids + 3] - 1111 * '0';
200 		i--;
201 		ids += 4;
202 	}
203 
204 	pbd = &d;
205 	powtwo = 0;
206 
207 	/* pre-scale to get the bits we want into the integer part */
208 	if (exp < 0) {
209 		/* i is a lower bound on log10(x) */
210 		i = exp + ndigs - 1;
211 		if (i <= 0 || ((i * 217705) >> 16) < sigbits + 2) {
212 			/*
213 			 * Scale by 2^(sigbits + 2 + u) where
214 			 * u is an upper bound on -log2(x).
215 			 */
216 			powtwo = sigbits + 2;
217 			if (i < 0)
218 				powtwo += ((-i * 217706) + 65535) >> 16;
219 			else if (i > 0)
220 				powtwo -= (i * 217705) >> 16;
221 			/*
222 			 * Take sigdigits large enough to get
223 			 * all integral digits correct.
224 			 */
225 			sigdigits = i + 1 + (((powtwo * 19729) + 65535) >> 16);
226 			__big_float_times_power(&d, 2, powtwo, sigdigits, &pbd);
227 		}
228 	}
229 
230 	/* convert to base 2^16 */
231 	b.bsize = _BIG_FLOAT_SIZE;
232 	pbb = &b;
233 	__big_decimal_to_big_binary(pbd, sigbits + 2, &pbb, &sticky);
234 
235 	/* adjust pbb->bexponent based on the scale factor above */
236 	pbb->bexponent -= powtwo;
237 
238 	/* convert to unpacked */
239 	ix = 0;
240 	for (i = pbb->blength - 1; i > 0 && ix < 5; i -= 2) {
241 		px->significand[ix++] = (pbb->bsignificand[i] << 16) |
242 		    pbb->bsignificand[i - 1];
243 	}
244 	if (ix < 5) {
245 		/* pad with zeroes */
246 		if (i == 0)
247 			px->significand[ix++] = pbb->bsignificand[i] << 16;
248 		while (ix < 5)
249 			px->significand[ix++] = 0;
250 	} else {
251 		/* truncate and set a sticky bit if necessary */
252 		while (i >= 0 && pbb->bsignificand[i] == 0)
253 			i--;
254 		if (i >= 0)
255 			px->significand[4] |= 1;
256 	}
257 	if (sticky | pd->more)
258 		px->significand[4] |= 1;
259 	px->exponent = pbb->bexponent + (pbb->blength << 4) - 1;
260 
261 	/* normalize so the most significant bit is set */
262 	while (px->significand[0] < 0x80000000u) {
263 		px->significand[0] = (px->significand[0] << 1) |
264 		    (px->significand[1] >> 31);
265 		px->significand[1] = (px->significand[1] << 1) |
266 		    (px->significand[2] >> 31);
267 		px->significand[2] = (px->significand[2] << 1) |
268 		    (px->significand[3] >> 31);
269 		px->significand[3] = (px->significand[3] << 1) |
270 		    (px->significand[4] >> 31);
271 		px->significand[4] <<= 1;
272 		px->exponent--;
273 	}
274 
275 	if (pbd != &d)
276 		(void) free((void *)pbd);
277 	if (pbb != &b)
278 		(void) free((void *)pbb);
279 }
280 
281 /*
282  * Convert a string s consisting of n <= 18 ASCII decimal digits
283  * to an integer value in double precision format, and set *pe
284  * to the number of rounding errors incurred (0 or 1).
285  */
286 static double
287 __digits_to_double(char *s, int n, int *pe)
288 {
289 	int	i, acc;
290 	double	t, th, tl;
291 
292 	if (n <= 9) {
293 		acc = s[0] - '0';
294 		for (i = 1; i < n; i++) {
295 			/* acc <- 10 * acc + next digit */
296 			acc = (acc << 1) + (acc << 3) + s[i] - '0';
297 		}
298 		t = (double)acc;
299 		*pe = 0;
300 	} else {
301 		acc = s[0] - '0';
302 		for (i = 1; i < (n - 9); i++) {
303 			/* acc <- 10 * acc + next digit */
304 			acc = (acc << 1) + (acc << 3) + s[i] - '0';
305 		}
306 		th = 1.0e9 * (double)acc;	/* this will be exact */
307 		acc = s[n - 9] - '0';
308 		for (i = n - 8; i < n; i++) {
309 			/* acc <- 10 * acc + next digit */
310 			acc = (acc << 1) + (acc << 3) + s[i] - '0';
311 		}
312 		tl = (double)acc;
313 		/* add and indicate whether or not the sum is exact */
314 		t = th + tl;
315 		*pe = ((t - th) == tl)? 0 : 1;
316 	}
317 	return (t);
318 }
319 
320 static union {
321 	int	i[2];
322 	double	d;
323 } C[] = {
324 #ifdef _LITTLE_ENDIAN
325 	{ 0x00000000, 0x3cc40000 },
326 #else
327 	{ 0x3cc40000, 0x00000000 },	/* 5 * 2^-53 */
328 #endif
329 };
330 
331 #define	five2m53	C[0].d
332 
333 static int
334 __fast_decimal_to_single(single *px, decimal_mode *pm, decimal_record *pd,
335     fp_exception_field_type *ps)
336 {
337 	double			dds, delta, ddsplus, ddsminus, df1;
338 	int			n, exp, rounded, e;
339 	float			f1, f2;
340 	__ieee_flags_type	fb;
341 
342 	if (pm->rd != fp_nearest)
343 		return (0);
344 
345 	exp = pd->exponent;
346 	if (pd->ndigits <= 18) {
347 		rounded = 0;
348 		n = pd->ndigits;
349 	} else {
350 		rounded = 1;
351 		n = 18;
352 		exp += pd->ndigits - 18;
353 	}
354 	/*
355 	 * exp must be in the range of the table, and the result
356 	 * must not underflow or overflow.
357 	 */
358 	if (exp < -__TBL_TENS_MAX || exp + n < -36 || exp + n > 38)
359 		return (0);
360 
361 	__get_ieee_flags(&fb);
362 	dds = __digits_to_double(pd->ds, n, &e);
363 	if (e != 0)
364 		rounded = 1;
365 	if (exp > 0) {
366 		/* small positive exponent */
367 		if (exp > __TBL_TENS_EXACT)
368 			rounded = 1;
369 		if (rounded) {
370 			dds *= __tbl_tens[exp];
371 		} else {
372 			dds = __mul_set(dds, __tbl_tens[exp], &e);
373 			if (e)
374 				rounded = 1;
375 		}
376 	} else if (exp < 0) {
377 		/* small negative exponent */
378 		if (-exp > __TBL_TENS_EXACT)
379 			rounded = 1;
380 		if (rounded) {
381 			dds /= __tbl_tens[-exp];
382 		} else {
383 			dds = __div_set(dds, __tbl_tens[-exp], &e);
384 			if (e)
385 				rounded = 1;
386 		}
387 	}
388 
389 	/*
390 	 * At this point dds may have four rounding errors due to
391 	 * (i) truncation of pd->ds to 18 digits, (ii) inexact con-
392 	 * version of pd->ds to binary, (iii) scaling by a power of
393 	 * ten that is not exactly representable, and (iv) roundoff
394 	 * error in the multiplication.  Below we will incur one more
395 	 * rounding error when we add or subtract delta and dds.  We
396 	 * construct delta so that even after this last rounding error,
397 	 * ddsplus is an upper bound on the exact value and ddsminus
398 	 * is a lower bound.  Then as long as both of these quantities
399 	 * round to the same single precision number, that number
400 	 * will be the correctly rounded single precision result.
401 	 * (If any rounding errors have been committed, then we must
402 	 * also be certain that the result can't be exact.)
403 	 */
404 	delta = five2m53 * dds;
405 	ddsplus = dds + delta;
406 	ddsminus = dds - delta;
407 	f1 = (float)(ddsplus);
408 	f2 = (float)(ddsminus);
409 	df1 = f1;
410 	__set_ieee_flags(&fb);
411 	if (f1 != f2)
412 		return (0);
413 	if (rounded) {
414 		/*
415 		 * If ddsminus <= f1 <= ddsplus, the result might be
416 		 * exact; we have to convert the long way to be sure.
417 		 */
418 		if (ddsminus <= df1 && df1 <= ddsplus)
419 			return (0);
420 		*ps = (1 << fp_inexact);
421 	} else {
422 		*ps = (df1 == dds)? 0 : (1 << fp_inexact);
423 	}
424 	*px = (pd->sign)? -f1 : f1;
425 	return (1);
426 }
427 
428 /*
429  * Attempt conversion to double using floating-point arithmetic.
430  * Return 1 if it works (at most one rounding error), 0 if it doesn't.
431  */
432 static int
433 __fast_decimal_to_double(double *px, decimal_mode *pm, decimal_record *pd,
434     fp_exception_field_type *ps)
435 {
436 	double			dds;
437 	int			e;
438 	__ieee_flags_type	fb;
439 
440 	if (pm->rd != fp_nearest || pd->ndigits > 18 || pd->exponent
441 	    > __TBL_TENS_EXACT || pd->exponent < -__TBL_TENS_EXACT)
442 		return (0);
443 
444 	__get_ieee_flags(&fb);
445 	dds = __digits_to_double(pd->ds, pd->ndigits, &e);
446 	if (e != 0) {
447 		__set_ieee_flags(&fb);
448 		return (0);
449 	}
450 	if (pd->exponent > 0)
451 		dds = __mul_set(dds, __tbl_tens[pd->exponent], &e);
452 	else if (pd->exponent < 0)
453 		dds = __div_set(dds, __tbl_tens[-pd->exponent], &e);
454 	*px = (pd->sign)? -dds : dds;
455 	*ps = (e)? (1 << fp_inexact) : 0;
456 	__set_ieee_flags(&fb);
457 	return (1);
458 }
459 
460 /* PUBLIC FUNCTIONS */
461 
462 /*
463  * The following routines convert the decimal record *pd to a floating
464  * point value *px observing the rounding mode specified in pm->rd and
465  * passing back any floating point exceptions that occur in *ps.
466  *
467  * pd->sign and pd->fpclass are always taken into account.  pd->exponent
468  * and pd->ds are used when pd->fpclass is fp_normal or fp_subnormal.
469  * In these cases, pd->ds must contain a null-terminated string of one
470  * or more ASCII digits, the first of which is not zero, and pd->ndigits
471  * must equal the length of this string.  If m is the integer represented
472  * by the string pd->ds, then *px will be set to a correctly rounded
473  * approximation to
474  *
475  *   -1**(pd->sign) * m * 10**(pd->exponent)
476  *
477  * (If pd->more != 0 then additional nonzero digits are assumed to follow
478  * those in pd->ds, so m is effectively replaced by m + epsilon in the
479  * expression above.)
480  *
481  * For example, if pd->exponent == -2 and pd->ds holds "1234", then *px
482  * will be a correctly rounded approximation to 12.34.
483  *
484  * Note that the only mode that matters is the rounding direction pm->rd;
485  * pm->df and pm->ndigits are never used.
486  */
487 
488 /* maximum decimal exponent we need to consider */
489 #define	SINGLE_MAXE	  47
490 #define	DOUBLE_MAXE	 326
491 #define	EXTENDED_MAXE	4968
492 #define	QUAD_MAXE	4968
493 
494 void
495 decimal_to_single(single *px, decimal_mode *pm, decimal_record *pd,
496     fp_exception_field_type *ps)
497 {
498 	single_equivalence	*kluge;
499 	unpacked		u;
500 	fp_exception_field_type	ef;
501 	int			i;
502 
503 	/* special values */
504 	kluge = (single_equivalence *)px;
505 	switch (pd->fpclass) {
506 	case fp_zero:
507 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
508 		kluge->f.msw.exponent = 0;
509 		kluge->f.msw.significand = 0;
510 		*ps = 0;
511 		return;
512 
513 	case fp_infinity:
514 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
515 		kluge->f.msw.exponent = 0xff;
516 		kluge->f.msw.significand = 0;
517 		*ps = 0;
518 		return;
519 
520 	case fp_quiet:
521 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
522 		kluge->f.msw.exponent = 0xff;
523 		kluge->f.msw.significand = 0x7fffff;
524 		*ps = 0;
525 		return;
526 
527 	case fp_signaling:
528 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
529 		kluge->f.msw.exponent = 0xff;
530 		kluge->f.msw.significand = 0x3fffff;
531 		*ps = 0;
532 		return;
533 	}
534 
535 	/* numeric values */
536 	ef = 0;
537 	if (pd->exponent + pd->ndigits > SINGLE_MAXE) {
538 		/* result must overflow */
539 		u.sign = (pd->sign != 0);
540 		u.fpclass = fp_normal;
541 		u.exponent = 0x000fffff;
542 		u.significand[0] = 0x80000000;
543 		for (i = 1; i < UNPACKED_SIZE; i++)
544 			u.significand[i] = 0;
545 	} else if (pd->exponent + pd->ndigits < -SINGLE_MAXE) {
546 		/* result must underflow completely */
547 		u.sign = (pd->sign != 0);
548 		u.fpclass = fp_normal;
549 		u.exponent = -0x000fffff;
550 		u.significand[0] = 0x80000000;
551 		for (i = 1; i < UNPACKED_SIZE; i++)
552 			u.significand[i] = 0;
553 	} else {
554 		/* result may be in range */
555 		if (__fast_decimal_to_single(px, pm, pd, &ef) == 1) {
556 			*ps = ef;
557 			if (ef != 0)
558 				__base_conversion_set_exception(ef);
559 			return;
560 		}
561 		__decimal_to_unpacked(&u, pd, 24);
562 	}
563 	__pack_single(&u, px, pm->rd, &ef);
564 	*ps = ef;
565 	if (ef != 0)
566 		__base_conversion_set_exception(ef);
567 }
568 
569 void
570 decimal_to_double(double *px, decimal_mode *pm, decimal_record *pd,
571     fp_exception_field_type *ps)
572 {
573 	double_equivalence	*kluge;
574 	unpacked		u;
575 	fp_exception_field_type	ef;
576 	int			i;
577 
578 	/* special values */
579 	kluge = (double_equivalence *)px;
580 	switch (pd->fpclass) {
581 	case fp_zero:
582 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
583 		kluge->f.msw.exponent = 0;
584 		kluge->f.msw.significand = 0;
585 		kluge->f.significand2 = 0;
586 		*ps = 0;
587 		return;
588 
589 	case fp_infinity:
590 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
591 		kluge->f.msw.exponent = 0x7ff;
592 		kluge->f.msw.significand = 0;
593 		kluge->f.significand2 = 0;
594 		*ps = 0;
595 		return;
596 
597 	case fp_quiet:
598 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
599 		kluge->f.msw.exponent = 0x7ff;
600 		kluge->f.msw.significand = 0xfffff;
601 		kluge->f.significand2 = 0xffffffff;
602 		*ps = 0;
603 		return;
604 
605 	case fp_signaling:
606 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
607 		kluge->f.msw.exponent = 0x7ff;
608 		kluge->f.msw.significand = 0x7ffff;
609 		kluge->f.significand2 = 0xffffffff;
610 		*ps = 0;
611 		return;
612 	}
613 
614 	/* numeric values */
615 	ef = 0;
616 	if (pd->exponent + pd->ndigits > DOUBLE_MAXE) {
617 		/* result must overflow */
618 		u.sign = (pd->sign != 0);
619 		u.fpclass = fp_normal;
620 		u.exponent = 0x000fffff;
621 		u.significand[0] = 0x80000000;
622 		for (i = 1; i < UNPACKED_SIZE; i++)
623 			u.significand[i] = 0;
624 	} else if (pd->exponent + pd->ndigits < -DOUBLE_MAXE) {
625 		/* result must underflow completely */
626 		u.sign = (pd->sign != 0);
627 		u.fpclass = fp_normal;
628 		u.exponent = -0x000fffff;
629 		u.significand[0] = 0x80000000;
630 		for (i = 1; i < UNPACKED_SIZE; i++)
631 			u.significand[i] = 0;
632 	} else {
633 		/* result may be in range */
634 		if (__fast_decimal_to_double(px, pm, pd, &ef) == 1) {
635 			*ps = ef;
636 			if (ef != 0)
637 				__base_conversion_set_exception(ef);
638 			return;
639 		}
640 		__decimal_to_unpacked(&u, pd, 53);
641 	}
642 	__pack_double(&u, px, pm->rd, &ef);
643 	*ps = ef;
644 	if (ef != 0)
645 		__base_conversion_set_exception(ef);
646 }
647 
648 void
649 decimal_to_extended(extended *px, decimal_mode *pm, decimal_record *pd,
650     fp_exception_field_type *ps)
651 {
652 	extended_equivalence	*kluge;
653 	unpacked		u;
654 	double_equivalence	dd;
655 	fp_exception_field_type ef;
656 	int			i;
657 
658 	/* special values */
659 	kluge = (extended_equivalence *)px;
660 	switch (pd->fpclass) {
661 	case fp_zero:
662 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
663 		kluge->f.msw.exponent = 0;
664 		kluge->f.significand = 0;
665 		kluge->f.significand2 = 0;
666 		*ps = 0;
667 		return;
668 
669 	case fp_infinity:
670 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
671 		kluge->f.msw.exponent = 0x7fff;
672 		kluge->f.significand = 0x80000000;
673 		kluge->f.significand2 = 0;
674 		*ps = 0;
675 		return;
676 
677 	case fp_quiet:
678 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
679 		kluge->f.msw.exponent = 0x7fff;
680 		kluge->f.significand = 0xffffffff;
681 		kluge->f.significand2 = 0xffffffff;
682 		*ps = 0;
683 		return;
684 
685 	case fp_signaling:
686 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
687 		kluge->f.msw.exponent = 0x7fff;
688 		kluge->f.significand = 0xbfffffff;
689 		kluge->f.significand2 = 0xffffffff;
690 		*ps = 0;
691 		return;
692 	}
693 
694 	/* numeric values */
695 	ef = 0;
696 	if (pd->exponent + pd->ndigits > EXTENDED_MAXE) {
697 		/* result must overflow */
698 		u.sign = (pd->sign != 0);
699 		u.fpclass = fp_normal;
700 		u.exponent = 0x000fffff;
701 		u.significand[0] = 0x80000000;
702 		for (i = 1; i < UNPACKED_SIZE; i++)
703 			u.significand[i] = 0;
704 	} else if (pd->exponent + pd->ndigits < -EXTENDED_MAXE) {
705 		/* result must underflow completely */
706 		u.sign = (pd->sign != 0);
707 		u.fpclass = fp_normal;
708 		u.exponent = -0x000fffff;
709 		u.significand[0] = 0x80000000;
710 		for (i = 1; i < UNPACKED_SIZE; i++)
711 			u.significand[i] = 0;
712 	} else {
713 		/* result may be in range */
714 		if (__fast_decimal_to_double(&dd.x, pm, pd, &ef) == 1 &&
715 		    ef == 0) {
716 			u.sign = dd.f.msw.sign;
717 			u.fpclass = fp_normal;
718 			u.exponent = dd.f.msw.exponent - DOUBLE_BIAS;
719 			u.significand[0] = ((0x100000 |
720 			    dd.f.msw.significand) << 11) |
721 			    (dd.f.significand2 >> 21);
722 			u.significand[1] = dd.f.significand2 << 11;
723 			for (i = 2; i < UNPACKED_SIZE; i++)
724 				u.significand[i] = 0;
725 		} else {
726 			__decimal_to_unpacked(&u, pd, 64);
727 		}
728 	}
729 	__pack_extended(&u, px, pm->rd, &ef);
730 	*ps = ef;
731 	if (ef != 0)
732 		__base_conversion_set_exception(ef);
733 }
734 
735 void
736 decimal_to_quadruple(quadruple *px, decimal_mode *pm, decimal_record *pd,
737     fp_exception_field_type *ps)
738 {
739 	quadruple_equivalence	*kluge;
740 	unpacked		u;
741 	double_equivalence	dd;
742 	fp_exception_field_type ef;
743 	int			i;
744 
745 	/* special values */
746 	kluge = (quadruple_equivalence *)px;
747 	switch (pd->fpclass) {
748 	case fp_zero:
749 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
750 		kluge->f.msw.exponent = 0;
751 		kluge->f.msw.significand = 0;
752 		kluge->f.significand2 = 0;
753 		kluge->f.significand3 = 0;
754 		kluge->f.significand4 = 0;
755 		*ps = 0;
756 		return;
757 
758 	case fp_infinity:
759 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
760 		kluge->f.msw.exponent = 0x7fff;
761 		kluge->f.msw.significand = 0;
762 		kluge->f.significand2 = 0;
763 		kluge->f.significand3 = 0;
764 		kluge->f.significand4 = 0;
765 		*ps = 0;
766 		return;
767 
768 	case fp_quiet:
769 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
770 		kluge->f.msw.exponent = 0x7fff;
771 		kluge->f.msw.significand = 0xffff;
772 		kluge->f.significand2 = 0xffffffff;
773 		kluge->f.significand3 = 0xffffffff;
774 		kluge->f.significand4 = 0xffffffff;
775 		*ps = 0;
776 		return;
777 
778 	case fp_signaling:
779 		kluge->f.msw.sign = (pd->sign)? 1 : 0;
780 		kluge->f.msw.exponent = 0x7fff;
781 		kluge->f.msw.significand = 0x7fff;
782 		kluge->f.significand2 = 0xffffffff;
783 		kluge->f.significand3 = 0xffffffff;
784 		kluge->f.significand4 = 0xffffffff;
785 		*ps = 0;
786 		return;
787 	}
788 
789 	/* numeric values */
790 	ef = 0;
791 	if (pd->exponent + pd->ndigits > QUAD_MAXE) {
792 		/* result must overflow */
793 		u.sign = (pd->sign != 0);
794 		u.fpclass = fp_normal;
795 		u.exponent = 0x000fffff;
796 		u.significand[0] = 0x80000000;
797 		for (i = 1; i < UNPACKED_SIZE; i++)
798 			u.significand[i] = 0;
799 	} else if (pd->exponent + pd->ndigits < -QUAD_MAXE) {
800 		/* result must underflow completely */
801 		u.sign = (pd->sign != 0);
802 		u.fpclass = fp_normal;
803 		u.exponent = -0x000fffff;
804 		u.significand[0] = 0x80000000;
805 		for (i = 1; i < UNPACKED_SIZE; i++)
806 			u.significand[i] = 0;
807 	} else {
808 		/* result may be in range */
809 		if (__fast_decimal_to_double(&dd.x, pm, pd, &ef) == 1 &&
810 		    ef == 0) {
811 			u.sign = dd.f.msw.sign;
812 			u.fpclass = fp_normal;
813 			u.exponent = dd.f.msw.exponent - DOUBLE_BIAS;
814 			u.significand[0] = ((0x100000 |
815 			    dd.f.msw.significand) << 11) |
816 			    (dd.f.significand2 >> 21);
817 			u.significand[1] = dd.f.significand2 << 11;
818 			for (i = 2; i < UNPACKED_SIZE; i++)
819 				u.significand[i] = 0;
820 		} else {
821 			__decimal_to_unpacked(&u, pd, 113);
822 		}
823 	}
824 	__pack_quadruple(&u, px, pm->rd, &ef);
825 	*ps = ef;
826 	if (ef != 0)
827 		__base_conversion_set_exception(ef);
828 }
829