sfmpy.c

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