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