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