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/dfmpy.c $Revision: 1.1 $
13 *
14 * Purpose:
15 * Double Precision Floating-point Multiply
16 *
17 * External Interfaces:
18 * dbl_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 "dbl_float.h"
31
32 /*
33 * Double Precision Floating-point Multiply
34 */
35
36 int
37 dbl_fmpy(
38 dbl_floating_point *srcptr1,
39 dbl_floating_point *srcptr2,
40 dbl_floating_point *dstptr,
41 unsigned int *status)
42 {
43 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
44 register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
45 register int dest_exponent, count;
46 register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
47 boolean is_tiny;
48
49 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
50 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
51
52 /*
53 * set sign bit of result
54 */
55 if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
56 Dbl_setnegativezerop1(resultp1);
57 else Dbl_setzerop1(resultp1);
58 /*
59 * check first operand for NaN's or infinity
60 */
61 if (Dbl_isinfinity_exponent(opnd1p1)) {
62 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
63 if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
64 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
65 /*
66 * invalid since operands are infinity
67 * and zero
68 */
69 if (Is_invalidtrap_enabled())
70 return(INVALIDEXCEPTION);
71 Set_invalidflag();
72 Dbl_makequietnan(resultp1,resultp2);
73 Dbl_copytoptr(resultp1,resultp2,dstptr);
74 return(NOEXCEPTION);
75 }
76 /*
77 * return infinity
78 */
79 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
80 Dbl_copytoptr(resultp1,resultp2,dstptr);
81 return(NOEXCEPTION);
82 }
83 }
84 else {
85 /*
86 * is NaN; signaling or quiet?
87 */
88 if (Dbl_isone_signaling(opnd1p1)) {
89 /* trap if INVALIDTRAP enabled */
90 if (Is_invalidtrap_enabled())
91 return(INVALIDEXCEPTION);
92 /* make NaN quiet */
93 Set_invalidflag();
94 Dbl_set_quiet(opnd1p1);
95 }
96 /*
97 * is second operand a signaling NaN?
98 */
99 else if (Dbl_is_signalingnan(opnd2p1)) {
100 /* trap if INVALIDTRAP enabled */
101 if (Is_invalidtrap_enabled())
102 return(INVALIDEXCEPTION);
103 /* make NaN quiet */
104 Set_invalidflag();
105 Dbl_set_quiet(opnd2p1);
106 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
107 return(NOEXCEPTION);
108 }
109 /*
110 * return quiet NaN
111 */
112 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
113 return(NOEXCEPTION);
114 }
115 }
116 /*
117 * check second operand for NaN's or infinity
118 */
119 if (Dbl_isinfinity_exponent(opnd2p1)) {
120 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
121 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
122 /* invalid since operands are zero & infinity */
123 if (Is_invalidtrap_enabled())
124 return(INVALIDEXCEPTION);
125 Set_invalidflag();
126 Dbl_makequietnan(opnd2p1,opnd2p2);
127 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
128 return(NOEXCEPTION);
129 }
130 /*
131 * return infinity
132 */
133 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
134 Dbl_copytoptr(resultp1,resultp2,dstptr);
135 return(NOEXCEPTION);
136 }
137 /*
138 * is NaN; signaling or quiet?
139 */
140 if (Dbl_isone_signaling(opnd2p1)) {
141 /* trap if INVALIDTRAP enabled */
142 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
143 /* make NaN quiet */
144 Set_invalidflag();
145 Dbl_set_quiet(opnd2p1);
146 }
147 /*
148 * return quiet NaN
149 */
150 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
151 return(NOEXCEPTION);
152 }
153 /*
154 * Generate exponent
155 */
156 dest_exponent = Dbl_exponent(opnd1p1) + Dbl_exponent(opnd2p1) -DBL_BIAS;
157
158 /*
159 * Generate mantissa
160 */
161 if (Dbl_isnotzero_exponent(opnd1p1)) {
162 /* set hidden bit */
163 Dbl_clear_signexponent_set_hidden(opnd1p1);
164 }
165 else {
166 /* check for zero */
167 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
168 Dbl_setzero_exponentmantissa(resultp1,resultp2);
169 Dbl_copytoptr(resultp1,resultp2,dstptr);
170 return(NOEXCEPTION);
171 }
172 /* is denormalized, adjust exponent */
173 Dbl_clear_signexponent(opnd1p1);
174 Dbl_leftshiftby1(opnd1p1,opnd1p2);
175 Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
176 }
177 /* opnd2 needs to have hidden bit set with msb in hidden bit */
178 if (Dbl_isnotzero_exponent(opnd2p1)) {
179 Dbl_clear_signexponent_set_hidden(opnd2p1);
180 }
181 else {
182 /* check for zero */
183 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
184 Dbl_setzero_exponentmantissa(resultp1,resultp2);
185 Dbl_copytoptr(resultp1,resultp2,dstptr);
186 return(NOEXCEPTION);
187 }
188 /* is denormalized; want to normalize */
189 Dbl_clear_signexponent(opnd2p1);
190 Dbl_leftshiftby1(opnd2p1,opnd2p2);
191 Dbl_normalize(opnd2p1,opnd2p2,dest_exponent);
192 }
193
194 /* Multiply two source mantissas together */
195
196 /* make room for guard bits */
197 Dbl_leftshiftby7(opnd2p1,opnd2p2);
198 Dbl_setzero(opnd3p1,opnd3p2);
199 /*
200 * Four bits at a time are inspected in each loop, and a
201 * simple shift and add multiply algorithm is used.
202 */
203 for (count=1;count<=DBL_P;count+=4) {
204 stickybit |= Dlow4p2(opnd3p2);
205 Dbl_rightshiftby4(opnd3p1,opnd3p2);
206 if (Dbit28p2(opnd1p2)) {
207 /* Twoword_add should be an ADDC followed by an ADD. */
208 Twoword_add(opnd3p1, opnd3p2, opnd2p1<<3 | opnd2p2>>29,
209 opnd2p2<<3);
210 }
211 if (Dbit29p2(opnd1p2)) {
212 Twoword_add(opnd3p1, opnd3p2, opnd2p1<<2 | opnd2p2>>30,
213 opnd2p2<<2);
214 }
215 if (Dbit30p2(opnd1p2)) {
216 Twoword_add(opnd3p1, opnd3p2, opnd2p1<<1 | opnd2p2>>31,
217 opnd2p2<<1);
218 }
219 if (Dbit31p2(opnd1p2)) {
220 Twoword_add(opnd3p1, opnd3p2, opnd2p1, opnd2p2);
221 }
222 Dbl_rightshiftby4(opnd1p1,opnd1p2);
223 }
224 if (Dbit3p1(opnd3p1)==0) {
225 Dbl_leftshiftby1(opnd3p1,opnd3p2);
226 }
227 else {
228 /* result mantissa >= 2. */
229 dest_exponent++;
230 }
231 /* check for denormalized result */
232 while (Dbit3p1(opnd3p1)==0) {
233 Dbl_leftshiftby1(opnd3p1,opnd3p2);
234 dest_exponent--;
235 }
236 /*
237 * check for guard, sticky and inexact bits
238 */
239 stickybit |= Dallp2(opnd3p2) << 25;
240 guardbit = (Dallp2(opnd3p2) << 24) >> 31;
241 inexact = guardbit | stickybit;
242
243 /* align result mantissa */
244 Dbl_rightshiftby8(opnd3p1,opnd3p2);
245
246 /*
247 * round result
248 */
249 if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) {
250 Dbl_clear_signexponent(opnd3p1);
251 switch (Rounding_mode()) {
252 case ROUNDPLUS:
253 if (Dbl_iszero_sign(resultp1))
254 Dbl_increment(opnd3p1,opnd3p2);
255 break;
256 case ROUNDMINUS:
257 if (Dbl_isone_sign(resultp1))
258 Dbl_increment(opnd3p1,opnd3p2);
259 break;
260 case ROUNDNEAREST:
261 if (guardbit) {
262 if (stickybit || Dbl_isone_lowmantissap2(opnd3p2))
263 Dbl_increment(opnd3p1,opnd3p2);
264 }
265 }
266 if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
267 }
268 Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);
269
270 /*
271 * Test for overflow
272 */
273 if (dest_exponent >= DBL_INFINITY_EXPONENT) {
274 /* trap if OVERFLOWTRAP enabled */
275 if (Is_overflowtrap_enabled()) {
276 /*
277 * Adjust bias of result
278 */
279 Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
280 Dbl_copytoptr(resultp1,resultp2,dstptr);
281 if (inexact)
282 if (Is_inexacttrap_enabled())
283 return (OVERFLOWEXCEPTION | INEXACTEXCEPTION);
284 else Set_inexactflag();
285 return (OVERFLOWEXCEPTION);
286 }
287 inexact = TRUE;
288 Set_overflowflag();
289 /* set result to infinity or largest number */
290 Dbl_setoverflow(resultp1,resultp2);
291 }
292 /*
293 * Test for underflow
294 */
295 else if (dest_exponent <= 0) {
296 /* trap if UNDERFLOWTRAP enabled */
297 if (Is_underflowtrap_enabled()) {
298 /*
299 * Adjust bias of result
300 */
301 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
302 Dbl_copytoptr(resultp1,resultp2,dstptr);
303 if (inexact)
304 if (Is_inexacttrap_enabled())
305 return (UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
306 else Set_inexactflag();
307 return (UNDERFLOWEXCEPTION);
308 }
309
310 /* Determine if should set underflow flag */
311 is_tiny = TRUE;
312 if (dest_exponent == 0 && inexact) {
313 switch (Rounding_mode()) {
314 case ROUNDPLUS:
315 if (Dbl_iszero_sign(resultp1)) {
316 Dbl_increment(opnd3p1,opnd3p2);
317 if (Dbl_isone_hiddenoverflow(opnd3p1))
318 is_tiny = FALSE;
319 Dbl_decrement(opnd3p1,opnd3p2);
320 }
321 break;
322 case ROUNDMINUS:
323 if (Dbl_isone_sign(resultp1)) {
324 Dbl_increment(opnd3p1,opnd3p2);
325 if (Dbl_isone_hiddenoverflow(opnd3p1))
326 is_tiny = FALSE;
327 Dbl_decrement(opnd3p1,opnd3p2);
328 }
329 break;
330 case ROUNDNEAREST:
331 if (guardbit && (stickybit ||
332 Dbl_isone_lowmantissap2(opnd3p2))) {
333 Dbl_increment(opnd3p1,opnd3p2);
334 if (Dbl_isone_hiddenoverflow(opnd3p1))
335 is_tiny = FALSE;
336 Dbl_decrement(opnd3p1,opnd3p2);
337 }
338 break;
339 }
340 }
341
342 /*
343 * denormalize result or set to signed zero
344 */
345 stickybit = inexact;
346 Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
347 stickybit,inexact);
348
349 /* return zero or smallest number */
350 if (inexact) {
351 switch (Rounding_mode()) {
352 case ROUNDPLUS:
353 if (Dbl_iszero_sign(resultp1)) {
354 Dbl_increment(opnd3p1,opnd3p2);
355 }
356 break;
357 case ROUNDMINUS:
358 if (Dbl_isone_sign(resultp1)) {
359 Dbl_increment(opnd3p1,opnd3p2);
360 }
361 break;
362 case ROUNDNEAREST:
363 if (guardbit && (stickybit ||
364 Dbl_isone_lowmantissap2(opnd3p2))) {
365 Dbl_increment(opnd3p1,opnd3p2);
366 }
367 break;
368 }
369 if (is_tiny) Set_underflowflag();
370 }
371 Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
372 }
373 else Dbl_set_exponent(resultp1,dest_exponent);
374 /* check for inexact */
375 Dbl_copytoptr(resultp1,resultp2,dstptr);
376 if (inexact) {
377 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
378 else Set_inexactflag();
379 }
380 return(NOEXCEPTION);
381 }