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