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