1 /*
2  * Copyright (c) 1992, 1993
3  *	The Regents of the University of California.  All rights reserved.
4  *
5  * This software was developed by the Computer Systems Engineering group
6  * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
7  * contributed to Berkeley.
8  *
9  * All advertising materials mentioning features or use of this software
10  * must display the following acknowledgement:
11  *	This product includes software developed by the University of
12  *	California, Lawrence Berkeley Laboratory.
13  *
14  * Redistribution and use in source and binary forms, with or without
15  * modification, are permitted provided that the following conditions
16  * are met:
17  * 1. Redistributions of source code must retain the above copyright
18  *    notice, this list of conditions and the following disclaimer.
19  * 2. Redistributions in binary form must reproduce the above copyright
20  *    notice, this list of conditions and the following disclaimer in the
21  *    documentation and/or other materials provided with the distribution.
22  * 3. Neither the name of the University nor the names of its contributors
23  *    may be used to endorse or promote products derived from this software
24  *    without specific prior written permission.
25  *
26  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
27  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
29  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
30  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
31  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
32  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
33  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
34  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
35  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
36  * SUCH DAMAGE.
37  *
38  *	@(#)fpu_implode.c	8.1 (Berkeley) 6/11/93
39  *	$NetBSD: fpu_implode.c,v 1.8 2001/08/26 05:44:46 eeh Exp $
40  */
41 
42 #include <sys/cdefs.h>
43 __FBSDID("$FreeBSD$");
44 
45 /*
46  * FPU subroutines: `implode' internal format numbers into the machine's
47  * `packed binary' format.
48  */
49 
50 #include <sys/param.h>
51 
52 #ifdef FPU_DEBUG
53 #include <stdio.h>
54 #endif
55 
56 #include <machine/frame.h>
57 #include <machine/fp.h>
58 #include <machine/fsr.h>
59 #include <machine/ieee.h>
60 #include <machine/instr.h>
61 
62 #include "fpu_arith.h"
63 #include "fpu_emu.h"
64 #include "fpu_extern.h"
65 #include "__sparc_utrap_private.h"
66 
67 static int fpround(struct fpemu *, struct fpn *);
68 static int toinf(struct fpemu *, int);
69 
70 /*
71  * Round a number (algorithm from Motorola MC68882 manual, modified for
72  * our internal format).  Set inexact exception if rounding is required.
73  * Return true iff we rounded up.
74  *
75  * After rounding, we discard the guard and round bits by shifting right
76  * 2 bits (a la fpu_shr(), but we do not bother with fp->fp_sticky).
77  * This saves effort later.
78  *
79  * Note that we may leave the value 2.0 in fp->fp_mant; it is the caller's
80  * responsibility to fix this if necessary.
81  */
82 static int
fpround(struct fpemu * fe,struct fpn * fp)83 fpround(struct fpemu *fe, struct fpn *fp)
84 {
85 	u_int m0, m1, m2, m3;
86 	int gr, s;
87 
88 	m0 = fp->fp_mant[0];
89 	m1 = fp->fp_mant[1];
90 	m2 = fp->fp_mant[2];
91 	m3 = fp->fp_mant[3];
92 	gr = m3 & 3;
93 	s = fp->fp_sticky;
94 
95 	/* mant >>= FP_NG */
96 	m3 = (m3 >> FP_NG) | (m2 << (32 - FP_NG));
97 	m2 = (m2 >> FP_NG) | (m1 << (32 - FP_NG));
98 	m1 = (m1 >> FP_NG) | (m0 << (32 - FP_NG));
99 	m0 >>= FP_NG;
100 
101 	if ((gr | s) == 0)	/* result is exact: no rounding needed */
102 		goto rounddown;
103 
104 	fe->fe_cx |= FSR_NX;	/* inexact */
105 
106 	/* Go to rounddown to round down; break to round up. */
107 	switch (FSR_GET_RD(fe->fe_fsr)) {
108 	case FSR_RD_N:
109 	default:
110 		/*
111 		 * Round only if guard is set (gr & 2).  If guard is set,
112 		 * but round & sticky both clear, then we want to round
113 		 * but have a tie, so round to even, i.e., add 1 iff odd.
114 		 */
115 		if ((gr & 2) == 0)
116 			goto rounddown;
117 		if ((gr & 1) || fp->fp_sticky || (m3 & 1))
118 			break;
119 		goto rounddown;
120 
121 	case FSR_RD_Z:
122 		/* Round towards zero, i.e., down. */
123 		goto rounddown;
124 
125 	case FSR_RD_NINF:
126 		/* Round towards -Inf: up if negative, down if positive. */
127 		if (fp->fp_sign)
128 			break;
129 		goto rounddown;
130 
131 	case FSR_RD_PINF:
132 		/* Round towards +Inf: up if positive, down otherwise. */
133 		if (!fp->fp_sign)
134 			break;
135 		goto rounddown;
136 	}
137 
138 	/* Bump low bit of mantissa, with carry. */
139 	FPU_ADDS(m3, m3, 1);
140 	FPU_ADDCS(m2, m2, 0);
141 	FPU_ADDCS(m1, m1, 0);
142 	FPU_ADDC(m0, m0, 0);
143 	fp->fp_mant[0] = m0;
144 	fp->fp_mant[1] = m1;
145 	fp->fp_mant[2] = m2;
146 	fp->fp_mant[3] = m3;
147 	return (1);
148 
149 rounddown:
150 	fp->fp_mant[0] = m0;
151 	fp->fp_mant[1] = m1;
152 	fp->fp_mant[2] = m2;
153 	fp->fp_mant[3] = m3;
154 	return (0);
155 }
156 
157 /*
158  * For overflow: return true if overflow is to go to +/-Inf, according
159  * to the sign of the overflowing result.  If false, overflow is to go
160  * to the largest magnitude value instead.
161  */
162 static int
toinf(struct fpemu * fe,int sign)163 toinf(struct fpemu *fe, int sign)
164 {
165 	int inf;
166 
167 	/* look at rounding direction */
168 	switch (FSR_GET_RD(fe->fe_fsr)) {
169 	default:
170 	case FSR_RD_N:		/* the nearest value is always Inf */
171 		inf = 1;
172 		break;
173 
174 	case FSR_RD_Z:		/* toward 0 => never towards Inf */
175 		inf = 0;
176 		break;
177 
178 	case FSR_RD_PINF:	/* toward +Inf iff positive */
179 		inf = sign == 0;
180 		break;
181 
182 	case FSR_RD_NINF:	/* toward -Inf iff negative */
183 		inf = sign;
184 		break;
185 	}
186 	return (inf);
187 }
188 
189 /*
190  * fpn -> int (int value returned as return value).
191  *
192  * N.B.: this conversion always rounds towards zero (this is a peculiarity
193  * of the SPARC instruction set).
194  */
195 u_int
__fpu_ftoi(fe,fp)196 __fpu_ftoi(fe, fp)
197 	struct fpemu *fe;
198 	struct fpn *fp;
199 {
200 	u_int i;
201 	int sign, exp;
202 
203 	sign = fp->fp_sign;
204 	switch (fp->fp_class) {
205 	case FPC_ZERO:
206 		return (0);
207 
208 	case FPC_NUM:
209 		/*
210 		 * If exp >= 2^32, overflow.  Otherwise shift value right
211 		 * into last mantissa word (this will not exceed 0xffffffff),
212 		 * shifting any guard and round bits out into the sticky
213 		 * bit.  Then ``round'' towards zero, i.e., just set an
214 		 * inexact exception if sticky is set (see round()).
215 		 * If the result is > 0x80000000, or is positive and equals
216 		 * 0x80000000, overflow; otherwise the last fraction word
217 		 * is the result.
218 		 */
219 		if ((exp = fp->fp_exp) >= 32)
220 			break;
221 		/* NB: the following includes exp < 0 cases */
222 		if (__fpu_shr(fp, FP_NMANT - 1 - exp) != 0)
223 			fe->fe_cx |= FSR_NX;
224 		i = fp->fp_mant[3];
225 		if (i >= ((u_int)0x80000000 + sign))
226 			break;
227 		return (sign ? -i : i);
228 
229 	default:		/* Inf, qNaN, sNaN */
230 		break;
231 	}
232 	/* overflow: replace any inexact exception with invalid */
233 	fe->fe_cx = (fe->fe_cx & ~FSR_NX) | FSR_NV;
234 	return (0x7fffffff + sign);
235 }
236 
237 /*
238  * fpn -> extended int (high bits of int value returned as return value).
239  *
240  * N.B.: this conversion always rounds towards zero (this is a peculiarity
241  * of the SPARC instruction set).
242  */
243 u_int
__fpu_ftox(fe,fp,res)244 __fpu_ftox(fe, fp, res)
245 	struct fpemu *fe;
246 	struct fpn *fp;
247 	u_int *res;
248 {
249 	u_int64_t i;
250 	int sign, exp;
251 
252 	sign = fp->fp_sign;
253 	switch (fp->fp_class) {
254 	case FPC_ZERO:
255 		i = 0;
256 		goto done;
257 
258 	case FPC_NUM:
259 		/*
260 		 * If exp >= 2^64, overflow.  Otherwise shift value
261 		 * right into last mantissa word (this will not exceed
262 		 * 0xffffffffffffffff), shifting any guard and round
263 		 * bits out into the sticky bit.  Then ``round'' towards
264 		 * zero, i.e., just set an inexact exception if sticky
265 		 * is set (see round()).
266 		 * If the result is > 0x8000000000000000, or is positive
267 		 * and equals 0x8000000000000000, overflow; otherwise
268 		 * the last fraction word is the result.
269 		 */
270 		if ((exp = fp->fp_exp) >= 64)
271 			break;
272 		/* NB: the following includes exp < 0 cases */
273 		if (__fpu_shr(fp, FP_NMANT - 1 - exp) != 0)
274 			fe->fe_cx |= FSR_NX;
275 		i = ((u_int64_t)fp->fp_mant[2]<<32)|fp->fp_mant[3];
276 		if (i >= ((u_int64_t)0x8000000000000000LL + sign))
277 			break;
278 		if (sign)
279 			i = -i;
280 		goto done;
281 
282 	default:		/* Inf, qNaN, sNaN */
283 		break;
284 	}
285 	/* overflow: replace any inexact exception with invalid */
286 	fe->fe_cx = (fe->fe_cx & ~FSR_NX) | FSR_NV;
287 	i = 0x7fffffffffffffffLL + sign;
288 done:
289 	res[1] = i & 0xffffffff;
290 	return (i >> 32);
291 }
292 
293 /*
294  * fpn -> single (32 bit single returned as return value).
295  * We assume <= 29 bits in a single-precision fraction (1.f part).
296  */
297 u_int
__fpu_ftos(fe,fp)298 __fpu_ftos(fe, fp)
299 	struct fpemu *fe;
300 	struct fpn *fp;
301 {
302 	u_int sign = fp->fp_sign << 31;
303 	int exp;
304 
305 #define	SNG_EXP(e)	((e) << SNG_FRACBITS)	/* makes e an exponent */
306 #define	SNG_MASK	(SNG_EXP(1) - 1)	/* mask for fraction */
307 
308 	/* Take care of non-numbers first. */
309 	if (ISNAN(fp)) {
310 		/*
311 		 * Preserve upper bits of NaN, per SPARC V8 appendix N.
312 		 * Note that fp->fp_mant[0] has the quiet bit set,
313 		 * even if it is classified as a signalling NaN.
314 		 */
315 		(void) __fpu_shr(fp, FP_NMANT - 1 - SNG_FRACBITS);
316 		exp = SNG_EXP_INFNAN;
317 		goto done;
318 	}
319 	if (ISINF(fp))
320 		return (sign | SNG_EXP(SNG_EXP_INFNAN));
321 	if (ISZERO(fp))
322 		return (sign);
323 
324 	/*
325 	 * Normals (including subnormals).  Drop all the fraction bits
326 	 * (including the explicit ``implied'' 1 bit) down into the
327 	 * single-precision range.  If the number is subnormal, move
328 	 * the ``implied'' 1 into the explicit range as well, and shift
329 	 * right to introduce leading zeroes.  Rounding then acts
330 	 * differently for normals and subnormals: the largest subnormal
331 	 * may round to the smallest normal (1.0 x 2^minexp), or may
332 	 * remain subnormal.  A number that is subnormal before rounding
333 	 * will signal an underflow if the result is inexact or if underflow
334 	 * traps are enabled.
335 	 *
336 	 * Rounding a normal, on the other hand, always produces another
337 	 * normal (although either way the result might be too big for
338 	 * single precision, and cause an overflow).  If rounding a
339 	 * normal produces 2.0 in the fraction, we need not adjust that
340 	 * fraction at all, since both 1.0 and 2.0 are zero under the
341 	 * fraction mask.
342 	 *
343 	 * Note that the guard and round bits vanish from the number after
344 	 * rounding.
345 	 */
346 	if ((exp = fp->fp_exp + SNG_EXP_BIAS) <= 0) {	/* subnormal */
347 		/* -NG for g,r; -SNG_FRACBITS-exp for fraction */
348 		(void) __fpu_shr(fp, FP_NMANT - FP_NG - SNG_FRACBITS - exp);
349 		if (fpround(fe, fp) && fp->fp_mant[3] == SNG_EXP(1)) {
350 			fe->fe_cx |= FSR_UF;
351 			return (sign | SNG_EXP(1) | 0);
352 		}
353 		if ((fe->fe_cx & FSR_NX) ||
354 		    (fe->fe_fsr & (FSR_UF << FSR_TEM_SHIFT)))
355 			fe->fe_cx |= FSR_UF;
356 		return (sign | SNG_EXP(0) | fp->fp_mant[3]);
357 	}
358 	/* -FP_NG for g,r; -1 for implied 1; -SNG_FRACBITS for fraction */
359 	(void) __fpu_shr(fp, FP_NMANT - FP_NG - 1 - SNG_FRACBITS);
360 #ifdef DIAGNOSTIC
361 	if ((fp->fp_mant[3] & SNG_EXP(1 << FP_NG)) == 0)
362 		__utrap_panic("fpu_ftos");
363 #endif
364 	if (fpround(fe, fp) && fp->fp_mant[3] == SNG_EXP(2))
365 		exp++;
366 	if (exp >= SNG_EXP_INFNAN) {
367 		/* overflow to inf or to max single */
368 		fe->fe_cx |= FSR_OF | FSR_NX;
369 		if (toinf(fe, sign))
370 			return (sign | SNG_EXP(SNG_EXP_INFNAN));
371 		return (sign | SNG_EXP(SNG_EXP_INFNAN - 1) | SNG_MASK);
372 	}
373 done:
374 	/* phew, made it */
375 	return (sign | SNG_EXP(exp) | (fp->fp_mant[3] & SNG_MASK));
376 }
377 
378 /*
379  * fpn -> double (32 bit high-order result returned; 32-bit low order result
380  * left in res[1]).  Assumes <= 61 bits in double precision fraction.
381  *
382  * This code mimics fpu_ftos; see it for comments.
383  */
384 u_int
__fpu_ftod(fe,fp,res)385 __fpu_ftod(fe, fp, res)
386 	struct fpemu *fe;
387 	struct fpn *fp;
388 	u_int *res;
389 {
390 	u_int sign = fp->fp_sign << 31;
391 	int exp;
392 
393 #define	DBL_EXP(e)	((e) << (DBL_FRACBITS & 31))
394 #define	DBL_MASK	(DBL_EXP(1) - 1)
395 
396 	if (ISNAN(fp)) {
397 		(void) __fpu_shr(fp, FP_NMANT - 1 - DBL_FRACBITS);
398 		exp = DBL_EXP_INFNAN;
399 		goto done;
400 	}
401 	if (ISINF(fp)) {
402 		sign |= DBL_EXP(DBL_EXP_INFNAN);
403 		goto zero;
404 	}
405 	if (ISZERO(fp)) {
406 zero:		res[1] = 0;
407 		return (sign);
408 	}
409 
410 	if ((exp = fp->fp_exp + DBL_EXP_BIAS) <= 0) {
411 		(void) __fpu_shr(fp, FP_NMANT - FP_NG - DBL_FRACBITS - exp);
412 		if (fpround(fe, fp) && fp->fp_mant[2] == DBL_EXP(1)) {
413 			fe->fe_cx |= FSR_UF;
414 			res[1] = 0;
415 			return (sign | DBL_EXP(1) | 0);
416 		}
417 		if ((fe->fe_cx & FSR_NX) ||
418 		    (fe->fe_fsr & (FSR_UF << FSR_TEM_SHIFT)))
419 			fe->fe_cx |= FSR_UF;
420 		exp = 0;
421 		goto done;
422 	}
423 	(void) __fpu_shr(fp, FP_NMANT - FP_NG - 1 - DBL_FRACBITS);
424 	if (fpround(fe, fp) && fp->fp_mant[2] == DBL_EXP(2))
425 		exp++;
426 	if (exp >= DBL_EXP_INFNAN) {
427 		fe->fe_cx |= FSR_OF | FSR_NX;
428 		if (toinf(fe, sign)) {
429 			res[1] = 0;
430 			return (sign | DBL_EXP(DBL_EXP_INFNAN) | 0);
431 		}
432 		res[1] = ~0;
433 		return (sign | DBL_EXP(DBL_EXP_INFNAN - 1) | DBL_MASK);
434 	}
435 done:
436 	res[1] = fp->fp_mant[3];
437 	return (sign | DBL_EXP(exp) | (fp->fp_mant[2] & DBL_MASK));
438 }
439 
440 /*
441  * fpn -> extended (32 bit high-order result returned; low-order fraction
442  * words left in res[1]..res[3]).  Like ftod, which is like ftos ... but
443  * our internal format *is* extended precision, plus 2 bits for guard/round,
444  * so we can avoid a small bit of work.
445  */
446 u_int
__fpu_ftoq(fe,fp,res)447 __fpu_ftoq(fe, fp, res)
448 	struct fpemu *fe;
449 	struct fpn *fp;
450 	u_int *res;
451 {
452 	u_int sign = fp->fp_sign << 31;
453 	int exp;
454 
455 #define	EXT_EXP(e)	((e) << (EXT_FRACBITS & 31))
456 #define	EXT_MASK	(EXT_EXP(1) - 1)
457 
458 	if (ISNAN(fp)) {
459 		(void) __fpu_shr(fp, 2);	/* since we are not rounding */
460 		exp = EXT_EXP_INFNAN;
461 		goto done;
462 	}
463 	if (ISINF(fp)) {
464 		sign |= EXT_EXP(EXT_EXP_INFNAN);
465 		goto zero;
466 	}
467 	if (ISZERO(fp)) {
468 zero:		res[1] = res[2] = res[3] = 0;
469 		return (sign);
470 	}
471 
472 	if ((exp = fp->fp_exp + EXT_EXP_BIAS) <= 0) {
473 		(void) __fpu_shr(fp, FP_NMANT - FP_NG - EXT_FRACBITS - exp);
474 		if (fpround(fe, fp) && fp->fp_mant[0] == EXT_EXP(1)) {
475 			fe->fe_cx |= FSR_UF;
476 			res[1] = res[2] = res[3] = 0;
477 			return (sign | EXT_EXP(1) | 0);
478 		}
479 		if ((fe->fe_cx & FSR_NX) ||
480 		    (fe->fe_fsr & (FSR_UF << FSR_TEM_SHIFT)))
481 			fe->fe_cx |= FSR_UF;
482 		exp = 0;
483 		goto done;
484 	}
485 	/* Since internal == extended, no need to shift here. */
486 	if (fpround(fe, fp) && fp->fp_mant[0] == EXT_EXP(2))
487 		exp++;
488 	if (exp >= EXT_EXP_INFNAN) {
489 		fe->fe_cx |= FSR_OF | FSR_NX;
490 		if (toinf(fe, sign)) {
491 			res[1] = res[2] = res[3] = 0;
492 			return (sign | EXT_EXP(EXT_EXP_INFNAN) | 0);
493 		}
494 		res[1] = res[2] = res[3] = ~0;
495 		return (sign | EXT_EXP(EXT_EXP_INFNAN - 1) | EXT_MASK);
496 	}
497 done:
498 	res[1] = fp->fp_mant[1];
499 	res[2] = fp->fp_mant[2];
500 	res[3] = fp->fp_mant[3];
501 	return (sign | EXT_EXP(exp) | (fp->fp_mant[0] & EXT_MASK));
502 }
503 
504 /*
505  * Implode an fpn, writing the result into the given space.
506  */
507 void
__fpu_implode(fe,fp,type,space)508 __fpu_implode(fe, fp, type, space)
509 	struct fpemu *fe;
510 	struct fpn *fp;
511 	int type;
512 	u_int *space;
513 {
514 
515 	switch (type) {
516 	case FTYPE_LNG:
517 		space[0] = __fpu_ftox(fe, fp, space);
518 		break;
519 
520 	case FTYPE_INT:
521 		space[0] = __fpu_ftoi(fe, fp);
522 		break;
523 
524 	case FTYPE_SNG:
525 		space[0] = __fpu_ftos(fe, fp);
526 		break;
527 
528 	case FTYPE_DBL:
529 		space[0] = __fpu_ftod(fe, fp, space);
530 		break;
531 
532 	case FTYPE_EXT:
533 		/* funky rounding precision options ?? */
534 		space[0] = __fpu_ftoq(fe, fp, space);
535 		break;
536 
537 	default:
538 		__utrap_panic("fpu_implode");
539 	}
540 	DPRINTF(FPE_REG, ("fpu_implode: %x %x %x %x\n",
541 		space[0], space[1], space[2], space[3]));
542 }
543