fpu_add.c revision 1.1 1 1.1 simonb /* $NetBSD: fpu_add.c,v 1.1 2001/06/13 06:01:46 simonb Exp $ */
2 1.1 simonb
3 1.1 simonb /*
4 1.1 simonb * Copyright (c) 1992, 1993
5 1.1 simonb * The Regents of the University of California. All rights reserved.
6 1.1 simonb *
7 1.1 simonb * This software was developed by the Computer Systems Engineering group
8 1.1 simonb * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
9 1.1 simonb * contributed to Berkeley.
10 1.1 simonb *
11 1.1 simonb * All advertising materials mentioning features or use of this software
12 1.1 simonb * must display the following acknowledgement:
13 1.1 simonb * This product includes software developed by the University of
14 1.1 simonb * California, Lawrence Berkeley Laboratory.
15 1.1 simonb *
16 1.1 simonb * Redistribution and use in source and binary forms, with or without
17 1.1 simonb * modification, are permitted provided that the following conditions
18 1.1 simonb * are met:
19 1.1 simonb * 1. Redistributions of source code must retain the above copyright
20 1.1 simonb * notice, this list of conditions and the following disclaimer.
21 1.1 simonb * 2. Redistributions in binary form must reproduce the above copyright
22 1.1 simonb * notice, this list of conditions and the following disclaimer in the
23 1.1 simonb * documentation and/or other materials provided with the distribution.
24 1.1 simonb * 3. All advertising materials mentioning features or use of this software
25 1.1 simonb * must display the following acknowledgement:
26 1.1 simonb * This product includes software developed by the University of
27 1.1 simonb * California, Berkeley and its contributors.
28 1.1 simonb * 4. Neither the name of the University nor the names of its contributors
29 1.1 simonb * may be used to endorse or promote products derived from this software
30 1.1 simonb * without specific prior written permission.
31 1.1 simonb *
32 1.1 simonb * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
33 1.1 simonb * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
34 1.1 simonb * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
35 1.1 simonb * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
36 1.1 simonb * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
37 1.1 simonb * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
38 1.1 simonb * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
39 1.1 simonb * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
40 1.1 simonb * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
41 1.1 simonb * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
42 1.1 simonb * SUCH DAMAGE.
43 1.1 simonb *
44 1.1 simonb * @(#)fpu_add.c 8.1 (Berkeley) 6/11/93
45 1.1 simonb */
46 1.1 simonb
47 1.1 simonb /*
48 1.1 simonb * Perform an FPU add (return x + y).
49 1.1 simonb *
50 1.1 simonb * To subtract, negate y and call add.
51 1.1 simonb */
52 1.1 simonb
53 1.1 simonb #include <sys/types.h>
54 1.1 simonb #if defined(DIAGNOSTIC)||defined(DEBUG)
55 1.1 simonb #include <sys/systm.h>
56 1.1 simonb #endif
57 1.1 simonb
58 1.1 simonb #include <machine/reg.h>
59 1.1 simonb #include <powerpc/instr.h>
60 1.1 simonb #include <machine/fpu.h>
61 1.1 simonb
62 1.1 simonb #include <powerpc/fpu/fpu_arith.h>
63 1.1 simonb #include <powerpc/fpu/fpu_emu.h>
64 1.1 simonb #include <powerpc/fpu/fpu_extern.h>
65 1.1 simonb
66 1.1 simonb struct fpn *
67 1.1 simonb fpu_add(struct fpemu *fe)
68 1.1 simonb {
69 1.1 simonb struct fpn *x = &fe->fe_f1, *y = &fe->fe_f2, *r;
70 1.1 simonb u_int r0, r1, r2, r3;
71 1.1 simonb int rd;
72 1.1 simonb
73 1.1 simonb /*
74 1.1 simonb * Put the `heavier' operand on the right (see fpu_emu.h).
75 1.1 simonb * Then we will have one of the following cases, taken in the
76 1.1 simonb * following order:
77 1.1 simonb *
78 1.1 simonb * - y = NaN. Implied: if only one is a signalling NaN, y is.
79 1.1 simonb * The result is y.
80 1.1 simonb * - y = Inf. Implied: x != NaN (is 0, number, or Inf: the NaN
81 1.1 simonb * case was taken care of earlier).
82 1.1 simonb * If x = -y, the result is NaN. Otherwise the result
83 1.1 simonb * is y (an Inf of whichever sign).
84 1.1 simonb * - y is 0. Implied: x = 0.
85 1.1 simonb * If x and y differ in sign (one positive, one negative),
86 1.1 simonb * the result is +0 except when rounding to -Inf. If same:
87 1.1 simonb * +0 + +0 = +0; -0 + -0 = -0.
88 1.1 simonb * - x is 0. Implied: y != 0.
89 1.1 simonb * Result is y.
90 1.1 simonb * - other. Implied: both x and y are numbers.
91 1.1 simonb * Do addition a la Hennessey & Patterson.
92 1.1 simonb */
93 1.1 simonb DPRINTF(FPE_REG, ("fpu_add:\n"));
94 1.1 simonb DUMPFPN(FPE_REG, x);
95 1.1 simonb DUMPFPN(FPE_REG, y);
96 1.1 simonb DPRINTF(FPE_REG, ("=>\n"));
97 1.1 simonb ORDER(x, y);
98 1.1 simonb if (ISNAN(y)) {
99 1.1 simonb fe->fe_cx |= FPSCR_VXSNAN;
100 1.1 simonb DUMPFPN(FPE_REG, y);
101 1.1 simonb return (y);
102 1.1 simonb }
103 1.1 simonb if (ISINF(y)) {
104 1.1 simonb if (ISINF(x) && x->fp_sign != y->fp_sign) {
105 1.1 simonb fe->fe_cx |= FPSCR_VXISI;
106 1.1 simonb return (fpu_newnan(fe));
107 1.1 simonb }
108 1.1 simonb DUMPFPN(FPE_REG, y);
109 1.1 simonb return (y);
110 1.1 simonb }
111 1.1 simonb rd = ((fe->fe_fpscr) & FPSCR_RN);
112 1.1 simonb if (ISZERO(y)) {
113 1.1 simonb if (rd != FSR_RD_RM) /* only -0 + -0 gives -0 */
114 1.1 simonb y->fp_sign &= x->fp_sign;
115 1.1 simonb else /* any -0 operand gives -0 */
116 1.1 simonb y->fp_sign |= x->fp_sign;
117 1.1 simonb DUMPFPN(FPE_REG, y);
118 1.1 simonb return (y);
119 1.1 simonb }
120 1.1 simonb if (ISZERO(x)) {
121 1.1 simonb DUMPFPN(FPE_REG, y);
122 1.1 simonb return (y);
123 1.1 simonb }
124 1.1 simonb /*
125 1.1 simonb * We really have two numbers to add, although their signs may
126 1.1 simonb * differ. Make the exponents match, by shifting the smaller
127 1.1 simonb * number right (e.g., 1.011 => 0.1011) and increasing its
128 1.1 simonb * exponent (2^3 => 2^4). Note that we do not alter the exponents
129 1.1 simonb * of x and y here.
130 1.1 simonb */
131 1.1 simonb r = &fe->fe_f3;
132 1.1 simonb r->fp_class = FPC_NUM;
133 1.1 simonb if (x->fp_exp == y->fp_exp) {
134 1.1 simonb r->fp_exp = x->fp_exp;
135 1.1 simonb r->fp_sticky = 0;
136 1.1 simonb } else {
137 1.1 simonb if (x->fp_exp < y->fp_exp) {
138 1.1 simonb /*
139 1.1 simonb * Try to avoid subtract case iii (see below).
140 1.1 simonb * This also guarantees that x->fp_sticky = 0.
141 1.1 simonb */
142 1.1 simonb SWAP(x, y);
143 1.1 simonb }
144 1.1 simonb /* now x->fp_exp > y->fp_exp */
145 1.1 simonb r->fp_exp = x->fp_exp;
146 1.1 simonb r->fp_sticky = fpu_shr(y, x->fp_exp - y->fp_exp);
147 1.1 simonb }
148 1.1 simonb r->fp_sign = x->fp_sign;
149 1.1 simonb if (x->fp_sign == y->fp_sign) {
150 1.1 simonb FPU_DECL_CARRY
151 1.1 simonb
152 1.1 simonb /*
153 1.1 simonb * The signs match, so we simply add the numbers. The result
154 1.1 simonb * may be `supernormal' (as big as 1.111...1 + 1.111...1, or
155 1.1 simonb * 11.111...0). If so, a single bit shift-right will fix it
156 1.1 simonb * (but remember to adjust the exponent).
157 1.1 simonb */
158 1.1 simonb /* r->fp_mant = x->fp_mant + y->fp_mant */
159 1.1 simonb FPU_ADDS(r->fp_mant[3], x->fp_mant[3], y->fp_mant[3]);
160 1.1 simonb FPU_ADDCS(r->fp_mant[2], x->fp_mant[2], y->fp_mant[2]);
161 1.1 simonb FPU_ADDCS(r->fp_mant[1], x->fp_mant[1], y->fp_mant[1]);
162 1.1 simonb FPU_ADDC(r0, x->fp_mant[0], y->fp_mant[0]);
163 1.1 simonb if ((r->fp_mant[0] = r0) >= FP_2) {
164 1.1 simonb (void) fpu_shr(r, 1);
165 1.1 simonb r->fp_exp++;
166 1.1 simonb }
167 1.1 simonb } else {
168 1.1 simonb FPU_DECL_CARRY
169 1.1 simonb
170 1.1 simonb /*
171 1.1 simonb * The signs differ, so things are rather more difficult.
172 1.1 simonb * H&P would have us negate the negative operand and add;
173 1.1 simonb * this is the same as subtracting the negative operand.
174 1.1 simonb * This is quite a headache. Instead, we will subtract
175 1.1 simonb * y from x, regardless of whether y itself is the negative
176 1.1 simonb * operand. When this is done one of three conditions will
177 1.1 simonb * hold, depending on the magnitudes of x and y:
178 1.1 simonb * case i) |x| > |y|. The result is just x - y,
179 1.1 simonb * with x's sign, but it may need to be normalized.
180 1.1 simonb * case ii) |x| = |y|. The result is 0 (maybe -0)
181 1.1 simonb * so must be fixed up.
182 1.1 simonb * case iii) |x| < |y|. We goofed; the result should
183 1.1 simonb * be (y - x), with the same sign as y.
184 1.1 simonb * We could compare |x| and |y| here and avoid case iii,
185 1.1 simonb * but that would take just as much work as the subtract.
186 1.1 simonb * We can tell case iii has occurred by an overflow.
187 1.1 simonb *
188 1.1 simonb * N.B.: since x->fp_exp >= y->fp_exp, x->fp_sticky = 0.
189 1.1 simonb */
190 1.1 simonb /* r->fp_mant = x->fp_mant - y->fp_mant */
191 1.1 simonb FPU_SET_CARRY(y->fp_sticky);
192 1.1 simonb FPU_SUBCS(r3, x->fp_mant[3], y->fp_mant[3]);
193 1.1 simonb FPU_SUBCS(r2, x->fp_mant[2], y->fp_mant[2]);
194 1.1 simonb FPU_SUBCS(r1, x->fp_mant[1], y->fp_mant[1]);
195 1.1 simonb FPU_SUBC(r0, x->fp_mant[0], y->fp_mant[0]);
196 1.1 simonb if (r0 < FP_2) {
197 1.1 simonb /* cases i and ii */
198 1.1 simonb if ((r0 | r1 | r2 | r3) == 0) {
199 1.1 simonb /* case ii */
200 1.1 simonb r->fp_class = FPC_ZERO;
201 1.1 simonb r->fp_sign = rd == FSR_RD_RM;
202 1.1 simonb return (r);
203 1.1 simonb }
204 1.1 simonb } else {
205 1.1 simonb /*
206 1.1 simonb * Oops, case iii. This can only occur when the
207 1.1 simonb * exponents were equal, in which case neither
208 1.1 simonb * x nor y have sticky bits set. Flip the sign
209 1.1 simonb * (to y's sign) and negate the result to get y - x.
210 1.1 simonb */
211 1.1 simonb #ifdef DIAGNOSTIC
212 1.1 simonb if (x->fp_exp != y->fp_exp || r->fp_sticky)
213 1.1 simonb panic("fpu_add");
214 1.1 simonb #endif
215 1.1 simonb r->fp_sign = y->fp_sign;
216 1.1 simonb FPU_SUBS(r3, 0, r3);
217 1.1 simonb FPU_SUBCS(r2, 0, r2);
218 1.1 simonb FPU_SUBCS(r1, 0, r1);
219 1.1 simonb FPU_SUBC(r0, 0, r0);
220 1.1 simonb }
221 1.1 simonb r->fp_mant[3] = r3;
222 1.1 simonb r->fp_mant[2] = r2;
223 1.1 simonb r->fp_mant[1] = r1;
224 1.1 simonb r->fp_mant[0] = r0;
225 1.1 simonb if (r0 < FP_1)
226 1.1 simonb fpu_norm(r);
227 1.1 simonb }
228 1.1 simonb DUMPFPN(FPE_REG, r);
229 1.1 simonb return (r);
230 1.1 simonb }
231