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