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