191174Stmm/*
291174Stmm * Copyright (c) 1992, 1993
391174Stmm *	The Regents of the University of California.  All rights reserved.
491174Stmm *
591174Stmm * This software was developed by the Computer Systems Engineering group
691174Stmm * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
791174Stmm * contributed to Berkeley.
891174Stmm *
991174Stmm * All advertising materials mentioning features or use of this software
1091174Stmm * must display the following acknowledgement:
1191174Stmm *	This product includes software developed by the University of
1291174Stmm *	California, Lawrence Berkeley Laboratory.
1391174Stmm *
1491174Stmm * Redistribution and use in source and binary forms, with or without
1591174Stmm * modification, are permitted provided that the following conditions
1691174Stmm * are met:
1791174Stmm * 1. Redistributions of source code must retain the above copyright
1891174Stmm *    notice, this list of conditions and the following disclaimer.
1991174Stmm * 2. Redistributions in binary form must reproduce the above copyright
2091174Stmm *    notice, this list of conditions and the following disclaimer in the
2191174Stmm *    documentation and/or other materials provided with the distribution.
2291174Stmm * 4. Neither the name of the University nor the names of its contributors
2391174Stmm *    may be used to endorse or promote products derived from this software
2491174Stmm *    without specific prior written permission.
2591174Stmm *
2691174Stmm * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
2791174Stmm * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
2891174Stmm * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
2991174Stmm * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
3091174Stmm * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
3191174Stmm * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
3291174Stmm * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
3391174Stmm * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
3491174Stmm * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
3591174Stmm * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
3691174Stmm * SUCH DAMAGE.
3791174Stmm *
3891174Stmm *	@(#)fpu_add.c	8.1 (Berkeley) 6/11/93
3992986Sobrien *	$NetBSD: fpu_add.c,v 1.3 1996/03/14 19:41:52 christos Exp $
4091174Stmm */
4191174Stmm
4292986Sobrien#include <sys/cdefs.h>
4392986Sobrien__FBSDID("$FreeBSD$");
4492986Sobrien
4591174Stmm/*
4691174Stmm * Perform an FPU add (return x + y).
4791174Stmm *
4891174Stmm * To subtract, negate y and call add.
4991174Stmm */
5091174Stmm
5191174Stmm#include <sys/param.h>
5291174Stmm
5391174Stmm#include <machine/frame.h>
5491174Stmm#include <machine/fp.h>
5591174Stmm#include <machine/fsr.h>
5691174Stmm#include <machine/instr.h>
5791174Stmm
5891174Stmm#include "fpu_arith.h"
5991174Stmm#include "fpu_emu.h"
6091174Stmm#include "fpu_extern.h"
6195587Sjake#include "__sparc_utrap_private.h"
6291174Stmm
6391174Stmmstruct fpn *
6491174Stmm__fpu_add(fe)
6592889Sobrien	struct fpemu *fe;
6691174Stmm{
6792889Sobrien	struct fpn *x = &fe->fe_f1, *y = &fe->fe_f2, *r;
6892889Sobrien	u_int r0, r1, r2, r3;
6992889Sobrien	int rd;
7091174Stmm
7191174Stmm	/*
7291174Stmm	 * Put the `heavier' operand on the right (see fpu_emu.h).
7391174Stmm	 * Then we will have one of the following cases, taken in the
7491174Stmm	 * following order:
7591174Stmm	 *
7691174Stmm	 *  - y = NaN.  Implied: if only one is a signalling NaN, y is.
7791174Stmm	 *	The result is y.
7891174Stmm	 *  - y = Inf.  Implied: x != NaN (is 0, number, or Inf: the NaN
7991174Stmm	 *    case was taken care of earlier).
8091174Stmm	 *	If x = -y, the result is NaN.  Otherwise the result
8191174Stmm	 *	is y (an Inf of whichever sign).
8291174Stmm	 *  - y is 0.  Implied: x = 0.
8391174Stmm	 *	If x and y differ in sign (one positive, one negative),
8491174Stmm	 *	the result is +0 except when rounding to -Inf.  If same:
8591174Stmm	 *	+0 + +0 = +0; -0 + -0 = -0.
8691174Stmm	 *  - x is 0.  Implied: y != 0.
8791174Stmm	 *	Result is y.
8891174Stmm	 *  - other.  Implied: both x and y are numbers.
8991174Stmm	 *	Do addition a la Hennessey & Patterson.
9091174Stmm	 */
9191174Stmm	ORDER(x, y);
9291174Stmm	if (ISNAN(y))
9391174Stmm		return (y);
9491174Stmm	if (ISINF(y)) {
9591174Stmm		if (ISINF(x) && x->fp_sign != y->fp_sign)
9691174Stmm			return (__fpu_newnan(fe));
9791174Stmm		return (y);
9891174Stmm	}
9991174Stmm	rd = FSR_GET_RD(fe->fe_fsr);
10091174Stmm	if (ISZERO(y)) {
10191174Stmm		if (rd != FSR_RD_NINF)	/* only -0 + -0 gives -0 */
10291174Stmm			y->fp_sign &= x->fp_sign;
10391174Stmm		else			/* any -0 operand gives -0 */
10491174Stmm			y->fp_sign |= x->fp_sign;
10591174Stmm		return (y);
10691174Stmm	}
10791174Stmm	if (ISZERO(x))
10891174Stmm		return (y);
10991174Stmm	/*
11091174Stmm	 * We really have two numbers to add, although their signs may
11191174Stmm	 * differ.  Make the exponents match, by shifting the smaller
11291174Stmm	 * number right (e.g., 1.011 => 0.1011) and increasing its
11391174Stmm	 * exponent (2^3 => 2^4).  Note that we do not alter the exponents
11491174Stmm	 * of x and y here.
11591174Stmm	 */
11691174Stmm	r = &fe->fe_f3;
11791174Stmm	r->fp_class = FPC_NUM;
11891174Stmm	if (x->fp_exp == y->fp_exp) {
11991174Stmm		r->fp_exp = x->fp_exp;
12091174Stmm		r->fp_sticky = 0;
12191174Stmm	} else {
12291174Stmm		if (x->fp_exp < y->fp_exp) {
12391174Stmm			/*
12491174Stmm			 * Try to avoid subtract case iii (see below).
12591174Stmm			 * This also guarantees that x->fp_sticky = 0.
12691174Stmm			 */
12791174Stmm			SWAP(x, y);
12891174Stmm		}
12991174Stmm		/* now x->fp_exp > y->fp_exp */
13091174Stmm		r->fp_exp = x->fp_exp;
13191174Stmm		r->fp_sticky = __fpu_shr(y, x->fp_exp - y->fp_exp);
13291174Stmm	}
13391174Stmm	r->fp_sign = x->fp_sign;
13491174Stmm	if (x->fp_sign == y->fp_sign) {
13591174Stmm		FPU_DECL_CARRY
13691174Stmm
13791174Stmm		/*
13891174Stmm		 * The signs match, so we simply add the numbers.  The result
13991174Stmm		 * may be `supernormal' (as big as 1.111...1 + 1.111...1, or
14091174Stmm		 * 11.111...0).  If so, a single bit shift-right will fix it
14191174Stmm		 * (but remember to adjust the exponent).
14291174Stmm		 */
14391174Stmm		/* r->fp_mant = x->fp_mant + y->fp_mant */
14491174Stmm		FPU_ADDS(r->fp_mant[3], x->fp_mant[3], y->fp_mant[3]);
14591174Stmm		FPU_ADDCS(r->fp_mant[2], x->fp_mant[2], y->fp_mant[2]);
14691174Stmm		FPU_ADDCS(r->fp_mant[1], x->fp_mant[1], y->fp_mant[1]);
14791174Stmm		FPU_ADDC(r0, x->fp_mant[0], y->fp_mant[0]);
14891174Stmm		if ((r->fp_mant[0] = r0) >= FP_2) {
14991174Stmm			(void) __fpu_shr(r, 1);
15091174Stmm			r->fp_exp++;
15191174Stmm		}
15291174Stmm	} else {
15391174Stmm		FPU_DECL_CARRY
15491174Stmm
15591174Stmm		/*
15691174Stmm		 * The signs differ, so things are rather more difficult.
15791174Stmm		 * H&P would have us negate the negative operand and add;
15891174Stmm		 * this is the same as subtracting the negative operand.
15991174Stmm		 * This is quite a headache.  Instead, we will subtract
16091174Stmm		 * y from x, regardless of whether y itself is the negative
16191174Stmm		 * operand.  When this is done one of three conditions will
16291174Stmm		 * hold, depending on the magnitudes of x and y:
16391174Stmm		 *   case i)   |x| > |y|.  The result is just x - y,
16491174Stmm		 *	with x's sign, but it may need to be normalized.
16591174Stmm		 *   case ii)  |x| = |y|.  The result is 0 (maybe -0)
16691174Stmm		 *	so must be fixed up.
16791174Stmm		 *   case iii) |x| < |y|.  We goofed; the result should
16891174Stmm		 *	be (y - x), with the same sign as y.
16991174Stmm		 * We could compare |x| and |y| here and avoid case iii,
17091174Stmm		 * but that would take just as much work as the subtract.
17191174Stmm		 * We can tell case iii has occurred by an overflow.
17291174Stmm		 *
17391174Stmm		 * N.B.: since x->fp_exp >= y->fp_exp, x->fp_sticky = 0.
17491174Stmm		 */
17591174Stmm		/* r->fp_mant = x->fp_mant - y->fp_mant */
17691174Stmm		FPU_SET_CARRY(y->fp_sticky);
17791174Stmm		FPU_SUBCS(r3, x->fp_mant[3], y->fp_mant[3]);
17891174Stmm		FPU_SUBCS(r2, x->fp_mant[2], y->fp_mant[2]);
17991174Stmm		FPU_SUBCS(r1, x->fp_mant[1], y->fp_mant[1]);
18091174Stmm		FPU_SUBC(r0, x->fp_mant[0], y->fp_mant[0]);
18191174Stmm		if (r0 < FP_2) {
18291174Stmm			/* cases i and ii */
18391174Stmm			if ((r0 | r1 | r2 | r3) == 0) {
18491174Stmm				/* case ii */
18591174Stmm				r->fp_class = FPC_ZERO;
18691174Stmm				r->fp_sign = rd == FSR_RD_NINF;
18791174Stmm				return (r);
18891174Stmm			}
18991174Stmm		} else {
19091174Stmm			/*
19191174Stmm			 * Oops, case iii.  This can only occur when the
19291174Stmm			 * exponents were equal, in which case neither
19391174Stmm			 * x nor y have sticky bits set.  Flip the sign
19491174Stmm			 * (to y's sign) and negate the result to get y - x.
19591174Stmm			 */
19691174Stmm#ifdef DIAGNOSTIC
19791174Stmm			if (x->fp_exp != y->fp_exp || r->fp_sticky)
19895587Sjake				__utrap_panic("fpu_add");
19991174Stmm#endif
20091174Stmm			r->fp_sign = y->fp_sign;
20191174Stmm			FPU_SUBS(r3, 0, r3);
20291174Stmm			FPU_SUBCS(r2, 0, r2);
20391174Stmm			FPU_SUBCS(r1, 0, r1);
20491174Stmm			FPU_SUBC(r0, 0, r0);
20591174Stmm		}
20691174Stmm		r->fp_mant[3] = r3;
20791174Stmm		r->fp_mant[2] = r2;
20891174Stmm		r->fp_mant[1] = r1;
20991174Stmm		r->fp_mant[0] = r0;
21091174Stmm		if (r0 < FP_1)
21191174Stmm			__fpu_norm(r);
21291174Stmm	}
21391174Stmm	return (r);
21491174Stmm}
215