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