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 * %sccs.include.redist.c%
15 *
16 * @(#)fpu_mul.c 8.1 (Berkeley) 06/11/93
17 *
18 * from: $Header: fpu_mul.c,v 1.3 92/11/26 01:39:50 torek Exp $
19 */
20
21 /*
22 * Perform an FPU multiply (return x * y).
23 */
24
25 #include <sys/types.h>
26
27 #include <machine/reg.h>
28
29 #include <sparc/fpu/fpu_arith.h>
30 #include <sparc/fpu/fpu_emu.h>
31
32 /*
33 * The multiplication algorithm for normal numbers is as follows:
34 *
35 * The fraction of the product is built in the usual stepwise fashion.
36 * Each step consists of shifting the accumulator right one bit
37 * (maintaining any guard bits) and, if the next bit in y is set,
38 * adding the multiplicand (x) to the accumulator. Then, in any case,
39 * we advance one bit leftward in y. Algorithmically:
40 *
41 * A = 0;
42 * for (bit = 0; bit < FP_NMANT; bit++) {
43 * sticky |= A & 1, A >>= 1;
44 * if (Y & (1 << bit))
45 * A += X;
46 * }
47 *
48 * (X and Y here represent the mantissas of x and y respectively.)
49 * The resultant accumulator (A) is the product's mantissa. It may
50 * be as large as 11.11111... in binary and hence may need to be
51 * shifted right, but at most one bit.
52 *
53 * Since we do not have efficient multiword arithmetic, we code the
54 * accumulator as four separate words, just like any other mantissa.
55 * We use local `register' variables in the hope that this is faster
56 * than memory. We keep x->fp_mant in locals for the same reason.
57 *
58 * In the algorithm above, the bits in y are inspected one at a time.
59 * We will pick them up 32 at a time and then deal with those 32, one
60 * at a time. Note, however, that we know several things about y:
61 *
62 * - the guard and round bits at the bottom are sure to be zero;
63 *
64 * - often many low bits are zero (y is often from a single or double
65 * precision source);
66 *
67 * - bit FP_NMANT-1 is set, and FP_1*2 fits in a word.
68 *
69 * We can also test for 32-zero-bits swiftly. In this case, the center
70 * part of the loop---setting sticky, shifting A, and not adding---will
71 * run 32 times without adding X to A. We can do a 32-bit shift faster
72 * by simply moving words. Since zeros are common, we optimize this case.
73 * Furthermore, since A is initially zero, we can omit the shift as well
74 * until we reach a nonzero word.
75 */
76 struct fpn *
fpu_mul(fe)77 fpu_mul(fe)
78 register struct fpemu *fe;
79 {
80 register struct fpn *x = &fe->fe_f1, *y = &fe->fe_f2;
81 register u_int a3, a2, a1, a0, x3, x2, x1, x0, bit, m;
82 register int sticky;
83 FPU_DECL_CARRY
84
85 /*
86 * Put the `heavier' operand on the right (see fpu_emu.h).
87 * Then we will have one of the following cases, taken in the
88 * following order:
89 *
90 * - y = NaN. Implied: if only one is a signalling NaN, y is.
91 * The result is y.
92 * - y = Inf. Implied: x != NaN (is 0, number, or Inf: the NaN
93 * case was taken care of earlier).
94 * If x = 0, the result is NaN. Otherwise the result
95 * is y, with its sign reversed if x is negative.
96 * - x = 0. Implied: y is 0 or number.
97 * The result is 0 (with XORed sign as usual).
98 * - other. Implied: both x and y are numbers.
99 * The result is x * y (XOR sign, multiply bits, add exponents).
100 */
101 ORDER(x, y);
102 if (ISNAN(y)) {
103 y->fp_sign ^= x->fp_sign;
104 return (y);
105 }
106 if (ISINF(y)) {
107 if (ISZERO(x))
108 return (fpu_newnan(fe));
109 y->fp_sign ^= x->fp_sign;
110 return (y);
111 }
112 if (ISZERO(x)) {
113 x->fp_sign ^= y->fp_sign;
114 return (x);
115 }
116
117 /*
118 * Setup. In the code below, the mask `m' will hold the current
119 * mantissa byte from y. The variable `bit' denotes the bit
120 * within m. We also define some macros to deal with everything.
121 */
122 x3 = x->fp_mant[3];
123 x2 = x->fp_mant[2];
124 x1 = x->fp_mant[1];
125 x0 = x->fp_mant[0];
126 sticky = a3 = a2 = a1 = a0 = 0;
127
128 #define ADD /* A += X */ \
129 FPU_ADDS(a3, a3, x3); \
130 FPU_ADDCS(a2, a2, x2); \
131 FPU_ADDCS(a1, a1, x1); \
132 FPU_ADDC(a0, a0, x0)
133
134 #define SHR1 /* A >>= 1, with sticky */ \
135 sticky |= a3 & 1, a3 = (a3 >> 1) | (a2 << 31), \
136 a2 = (a2 >> 1) | (a1 << 31), a1 = (a1 >> 1) | (a0 << 31), a0 >>= 1
137
138 #define SHR32 /* A >>= 32, with sticky */ \
139 sticky |= a3, a3 = a2, a2 = a1, a1 = a0, a0 = 0
140
141 #define STEP /* each 1-bit step of the multiplication */ \
142 SHR1; if (bit & m) { ADD; }; bit <<= 1
143
144 /*
145 * We are ready to begin. The multiply loop runs once for each
146 * of the four 32-bit words. Some words, however, are special.
147 * As noted above, the low order bits of Y are often zero. Even
148 * if not, the first loop can certainly skip the guard bits.
149 * The last word of y has its highest 1-bit in position FP_NMANT-1,
150 * so we stop the loop when we move past that bit.
151 */
152 if ((m = y->fp_mant[3]) == 0) {
153 /* SHR32; */ /* unneeded since A==0 */
154 } else {
155 bit = 1 << FP_NG;
156 do {
157 STEP;
158 } while (bit != 0);
159 }
160 if ((m = y->fp_mant[2]) == 0) {
161 SHR32;
162 } else {
163 bit = 1;
164 do {
165 STEP;
166 } while (bit != 0);
167 }
168 if ((m = y->fp_mant[1]) == 0) {
169 SHR32;
170 } else {
171 bit = 1;
172 do {
173 STEP;
174 } while (bit != 0);
175 }
176 m = y->fp_mant[0]; /* definitely != 0 */
177 bit = 1;
178 do {
179 STEP;
180 } while (bit <= m);
181
182 /*
183 * Done with mantissa calculation. Get exponent and handle
184 * 11.111...1 case, then put result in place. We reuse x since
185 * it already has the right class (FP_NUM).
186 */
187 m = x->fp_exp + y->fp_exp;
188 if (a0 >= FP_2) {
189 SHR1;
190 m++;
191 }
192 x->fp_sign ^= y->fp_sign;
193 x->fp_exp = m;
194 x->fp_sticky = sticky;
195 x->fp_mant[3] = a3;
196 x->fp_mant[2] = a2;
197 x->fp_mant[1] = a1;
198 x->fp_mant[0] = a0;
199 return (x);
200 }
201