1 /* $OpenBSD: sfmpy.c,v 1.5 2002/05/07 22:19:30 mickey Exp $ */
2 /*
3 (c) Copyright 1986 HEWLETT-PACKARD COMPANY
4 To anyone who acknowledges that this file is provided "AS IS"
5 without any express or implied warranty:
6 permission to use, copy, modify, and distribute this file
7 for any purpose is hereby granted without fee, provided that
8 the above copyright notice and this notice appears in all
9 copies, and that the name of Hewlett-Packard Company not be
10 used in advertising or publicity pertaining to distribution
11 of the software without specific, written prior permission.
12 Hewlett-Packard Company makes no representations about the
13 suitability of this software for any purpose.
14 */
15 /* @(#)sfmpy.c: Revision: 2.9.88.1 Date: 93/12/07 15:07:07 */
16
17 #include "float.h"
18 #include "sgl_float.h"
19
20 /*
21 * Single Precision Floating-point Multiply
22 */
23 int
sgl_fmpy(srcptr1,srcptr2,dstptr,status)24 sgl_fmpy(srcptr1,srcptr2,dstptr,status)
25 sgl_floating_point *srcptr1, *srcptr2, *dstptr;
26 unsigned int *status;
27 {
28 register unsigned int opnd1, opnd2, opnd3, result;
29 register int dest_exponent, count;
30 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
31 int is_tiny;
32
33 opnd1 = *srcptr1;
34 opnd2 = *srcptr2;
35 /*
36 * set sign bit of result
37 */
38 if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);
39 else Sgl_setzero(result);
40 /*
41 * check first operand for NaN's or infinity
42 */
43 if (Sgl_isinfinity_exponent(opnd1)) {
44 if (Sgl_iszero_mantissa(opnd1)) {
45 if (Sgl_isnotnan(opnd2)) {
46 if (Sgl_iszero_exponentmantissa(opnd2)) {
47 /*
48 * invalid since operands are infinity
49 * and zero
50 */
51 if (Is_invalidtrap_enabled())
52 return(INVALIDEXCEPTION);
53 Set_invalidflag();
54 Sgl_makequietnan(result);
55 *dstptr = result;
56 return(NOEXCEPTION);
57 }
58 /*
59 * return infinity
60 */
61 Sgl_setinfinity_exponentmantissa(result);
62 *dstptr = result;
63 return(NOEXCEPTION);
64 }
65 }
66 else {
67 /*
68 * is NaN; signaling or quiet?
69 */
70 if (Sgl_isone_signaling(opnd1)) {
71 /* trap if INVALIDTRAP enabled */
72 if (Is_invalidtrap_enabled())
73 return(INVALIDEXCEPTION);
74 /* make NaN quiet */
75 Set_invalidflag();
76 Sgl_set_quiet(opnd1);
77 }
78 /*
79 * is second operand a signaling NaN?
80 */
81 else if (Sgl_is_signalingnan(opnd2)) {
82 /* trap if INVALIDTRAP enabled */
83 if (Is_invalidtrap_enabled())
84 return(INVALIDEXCEPTION);
85 /* make NaN quiet */
86 Set_invalidflag();
87 Sgl_set_quiet(opnd2);
88 *dstptr = opnd2;
89 return(NOEXCEPTION);
90 }
91 /*
92 * return quiet NaN
93 */
94 *dstptr = opnd1;
95 return(NOEXCEPTION);
96 }
97 }
98 /*
99 * check second operand for NaN's or infinity
100 */
101 if (Sgl_isinfinity_exponent(opnd2)) {
102 if (Sgl_iszero_mantissa(opnd2)) {
103 if (Sgl_iszero_exponentmantissa(opnd1)) {
104 /* invalid since operands are zero & infinity */
105 if (Is_invalidtrap_enabled())
106 return(INVALIDEXCEPTION);
107 Set_invalidflag();
108 Sgl_makequietnan(opnd2);
109 *dstptr = opnd2;
110 return(NOEXCEPTION);
111 }
112 /*
113 * return infinity
114 */
115 Sgl_setinfinity_exponentmantissa(result);
116 *dstptr = result;
117 return(NOEXCEPTION);
118 }
119 /*
120 * is NaN; signaling or quiet?
121 */
122 if (Sgl_isone_signaling(opnd2)) {
123 /* trap if INVALIDTRAP enabled */
124 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
125
126 /* make NaN quiet */
127 Set_invalidflag();
128 Sgl_set_quiet(opnd2);
129 }
130 /*
131 * return quiet NaN
132 */
133 *dstptr = opnd2;
134 return(NOEXCEPTION);
135 }
136 /*
137 * Generate exponent
138 */
139 dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS;
140
141 /*
142 * Generate mantissa
143 */
144 if (Sgl_isnotzero_exponent(opnd1)) {
145 /* set hidden bit */
146 Sgl_clear_signexponent_set_hidden(opnd1);
147 }
148 else {
149 /* check for zero */
150 if (Sgl_iszero_mantissa(opnd1)) {
151 Sgl_setzero_exponentmantissa(result);
152 *dstptr = result;
153 return(NOEXCEPTION);
154 }
155 /* is denormalized, adjust exponent */
156 Sgl_clear_signexponent(opnd1);
157 Sgl_leftshiftby1(opnd1);
158 Sgl_normalize(opnd1,dest_exponent);
159 }
160 /* opnd2 needs to have hidden bit set with msb in hidden bit */
161 if (Sgl_isnotzero_exponent(opnd2)) {
162 Sgl_clear_signexponent_set_hidden(opnd2);
163 }
164 else {
165 /* check for zero */
166 if (Sgl_iszero_mantissa(opnd2)) {
167 Sgl_setzero_exponentmantissa(result);
168 *dstptr = result;
169 return(NOEXCEPTION);
170 }
171 /* is denormalized; want to normalize */
172 Sgl_clear_signexponent(opnd2);
173 Sgl_leftshiftby1(opnd2);
174 Sgl_normalize(opnd2,dest_exponent);
175 }
176
177 /* Multiply two source mantissas together */
178
179 Sgl_leftshiftby4(opnd2); /* make room for guard bits */
180 Sgl_setzero(opnd3);
181 /*
182 * Four bits at a time are inspected in each loop, and a
183 * simple shift and add multiply algorithm is used.
184 */
185 for (count=1;count<SGL_P;count+=4) {
186 stickybit |= Slow4(opnd3);
187 Sgl_rightshiftby4(opnd3);
188 if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3);
189 if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2);
190 if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1);
191 if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2);
192 Sgl_rightshiftby4(opnd1);
193 }
194 /* make sure result is left-justified */
195 if (Sgl_iszero_sign(opnd3)) {
196 Sgl_leftshiftby1(opnd3);
197 }
198 else {
199 /* result mantissa >= 2. */
200 dest_exponent++;
201 }
202 /* check for denormalized result */
203 while (Sgl_iszero_sign(opnd3)) {
204 Sgl_leftshiftby1(opnd3);
205 dest_exponent--;
206 }
207 /*
208 * check for guard, sticky and inexact bits
209 */
210 stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1);
211 guardbit = Sbit24(opnd3);
212 inexact = guardbit | stickybit;
213
214 /* re-align mantissa */
215 Sgl_rightshiftby8(opnd3);
216
217 /*
218 * round result
219 */
220 if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) {
221 Sgl_clear_signexponent(opnd3);
222 switch (Rounding_mode()) {
223 case ROUNDPLUS:
224 if (Sgl_iszero_sign(result))
225 Sgl_increment(opnd3);
226 break;
227 case ROUNDMINUS:
228 if (Sgl_isone_sign(result))
229 Sgl_increment(opnd3);
230 break;
231 case ROUNDNEAREST:
232 if (guardbit &&
233 (stickybit || Sgl_isone_lowmantissa(opnd3)))
234 Sgl_increment(opnd3);
235 break;
236 }
237 if (Sgl_isone_hidden(opnd3)) dest_exponent++;
238 }
239 Sgl_set_mantissa(result,opnd3);
240
241 /*
242 * Test for overflow
243 */
244 if (dest_exponent >= SGL_INFINITY_EXPONENT) {
245 /* trap if OVERFLOWTRAP enabled */
246 if (Is_overflowtrap_enabled()) {
247 /*
248 * Adjust bias of result
249 */
250 Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
251 *dstptr = result;
252 if (inexact) {
253 if (Is_inexacttrap_enabled())
254 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
255 else Set_inexactflag();
256 }
257 return(OVERFLOWEXCEPTION);
258 }
259 inexact = TRUE;
260 Set_overflowflag();
261 /* set result to infinity or largest number */
262 Sgl_setoverflow(result);
263 }
264 /*
265 * Test for underflow
266 */
267 else if (dest_exponent <= 0) {
268 /* trap if UNDERFLOWTRAP enabled */
269 if (Is_underflowtrap_enabled()) {
270 /*
271 * Adjust bias of result
272 */
273 Sgl_setwrapped_exponent(result,dest_exponent,unfl);
274 *dstptr = result;
275 if (inexact) {
276 if (Is_inexacttrap_enabled())
277 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
278 else Set_inexactflag();
279 }
280 return(UNDERFLOWEXCEPTION);
281 }
282
283 /* Determine if should set underflow flag */
284 is_tiny = TRUE;
285 if (dest_exponent == 0 && inexact) {
286 switch (Rounding_mode()) {
287 case ROUNDPLUS:
288 if (Sgl_iszero_sign(result)) {
289 Sgl_increment(opnd3);
290 if (Sgl_isone_hiddenoverflow(opnd3))
291 is_tiny = FALSE;
292 Sgl_decrement(opnd3);
293 }
294 break;
295 case ROUNDMINUS:
296 if (Sgl_isone_sign(result)) {
297 Sgl_increment(opnd3);
298 if (Sgl_isone_hiddenoverflow(opnd3))
299 is_tiny = FALSE;
300 Sgl_decrement(opnd3);
301 }
302 break;
303 case ROUNDNEAREST:
304 if (guardbit && (stickybit ||
305 Sgl_isone_lowmantissa(opnd3))) {
306 Sgl_increment(opnd3);
307 if (Sgl_isone_hiddenoverflow(opnd3))
308 is_tiny = FALSE;
309 Sgl_decrement(opnd3);
310 }
311 break;
312 }
313 }
314
315 /*
316 * denormalize result or set to signed zero
317 */
318 stickybit = inexact;
319 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
320
321 /* return zero or smallest number */
322 if (inexact) {
323 switch (Rounding_mode()) {
324 case ROUNDPLUS:
325 if (Sgl_iszero_sign(result))
326 Sgl_increment(opnd3);
327 break;
328 case ROUNDMINUS:
329 if (Sgl_isone_sign(result))
330 Sgl_increment(opnd3);
331 break;
332 case ROUNDNEAREST:
333 if (guardbit && (stickybit ||
334 Sgl_isone_lowmantissa(opnd3)))
335 Sgl_increment(opnd3);
336 break;
337 }
338 if (is_tiny) Set_underflowflag();
339 }
340 Sgl_set_exponentmantissa(result,opnd3);
341 }
342 else Sgl_set_exponent(result,dest_exponent);
343 *dstptr = result;
344
345 /* check for inexact */
346 if (inexact) {
347 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
348 else Set_inexactflag();
349 }
350 return(NOEXCEPTION);
351 }
352