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