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