xref: /netbsd/sys/arch/hppa/spmath/sfdiv.c (revision c4a72b64) 
1 /*	$NetBSD: sfdiv.c,v 1.1 2002/06/05 01:04:26 fredette Exp $	*/
2 
3 /*	$OpenBSD: sfdiv.c,v 1.4 2001/03/29 03:58:19 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 "../spmath/float.h"
45 #include "../spmath/sgl_float.h"
46 
47 /*
48  *  Single Precision Floating-point Divide
49  */
50 int
51 sgl_fdiv(srcptr1,srcptr2,dstptr,status)
52 
53 sgl_floating_point *srcptr1, *srcptr2, *dstptr;
54 unsigned int *status;
55 {
56 	register unsigned int opnd1, opnd2, opnd3, result;
57 	register int dest_exponent, count;
58 	register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
59 	int 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_isinfinity(opnd2)) {
75 					/*
76 					 * invalid since both operands
77 					 * are infinity
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 			/*
132 			 * return zero
133 			 */
134 			Sgl_setzero_exponentmantissa(result);
135 			*dstptr = result;
136 			return(NOEXCEPTION);
137 		}
138 		/*
139 		 * is NaN; signaling or quiet?
140 		 */
141 		if (Sgl_isone_signaling(opnd2)) {
142 			/* trap if INVALIDTRAP enabled */
143 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
144 			/* make NaN quiet */
145 			Set_invalidflag();
146 			Sgl_set_quiet(opnd2);
147 		}
148 		/*
149 		 * return quiet NaN
150 		 */
151 		*dstptr = opnd2;
152 		return(NOEXCEPTION);
153 	}
154 	/*
155 	 * check for division by zero
156 	 */
157 	if (Sgl_iszero_exponentmantissa(opnd2)) {
158 		if (Sgl_iszero_exponentmantissa(opnd1)) {
159 			/* invalid since both operands are zero */
160 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
161 			Set_invalidflag();
162 			Sgl_makequietnan(result);
163 			*dstptr = result;
164 			return(NOEXCEPTION);
165 		}
166 		if (Is_divisionbyzerotrap_enabled())
167 			return(DIVISIONBYZEROEXCEPTION);
168 		Set_divisionbyzeroflag();
169 		Sgl_setinfinity_exponentmantissa(result);
170 		*dstptr = result;
171 		return(NOEXCEPTION);
172 	}
173 	/*
174 	 * Generate exponent
175 	 */
176 	dest_exponent = Sgl_exponent(opnd1) - Sgl_exponent(opnd2) + SGL_BIAS;
177 
178 	/*
179 	 * Generate mantissa
180 	 */
181 	if (Sgl_isnotzero_exponent(opnd1)) {
182 		/* set hidden bit */
183 		Sgl_clear_signexponent_set_hidden(opnd1);
184 	}
185 	else {
186 		/* check for zero */
187 		if (Sgl_iszero_mantissa(opnd1)) {
188 			Sgl_setzero_exponentmantissa(result);
189 			*dstptr = result;
190 			return(NOEXCEPTION);
191 		}
192 		/* is denormalized; want to normalize */
193 		Sgl_clear_signexponent(opnd1);
194 		Sgl_leftshiftby1(opnd1);
195 		Sgl_normalize(opnd1,dest_exponent);
196 	}
197 	/* opnd2 needs to have hidden bit set with msb in hidden bit */
198 	if (Sgl_isnotzero_exponent(opnd2)) {
199 		Sgl_clear_signexponent_set_hidden(opnd2);
200 	}
201 	else {
202 		/* is denormalized; want to normalize */
203 		Sgl_clear_signexponent(opnd2);
204 		Sgl_leftshiftby1(opnd2);
205 		while(Sgl_iszero_hiddenhigh7mantissa(opnd2)) {
206 			Sgl_leftshiftby8(opnd2);
207 			dest_exponent += 8;
208 		}
209 		if(Sgl_iszero_hiddenhigh3mantissa(opnd2)) {
210 			Sgl_leftshiftby4(opnd2);
211 			dest_exponent += 4;
212 		}
213 		while(Sgl_iszero_hidden(opnd2)) {
214 			Sgl_leftshiftby1(opnd2);
215 			dest_exponent += 1;
216 		}
217 	}
218 
219 	/* Divide the source mantissas */
220 
221 	/*
222 	 * A non_restoring divide algorithm is used.
223 	 */
224 	Sgl_subtract(opnd1,opnd2,opnd1);
225 	Sgl_setzero(opnd3);
226 	for (count=1;count<=SGL_P && Sgl_all(opnd1);count++) {
227 		Sgl_leftshiftby1(opnd1);
228 		Sgl_leftshiftby1(opnd3);
229 		if (Sgl_iszero_sign(opnd1)) {
230 			Sgl_setone_lowmantissa(opnd3);
231 			Sgl_subtract(opnd1,opnd2,opnd1);
232 		}
233 		else Sgl_addition(opnd1,opnd2,opnd1);
234 	}
235 	if (count <= SGL_P) {
236 		Sgl_leftshiftby1(opnd3);
237 		Sgl_setone_lowmantissa(opnd3);
238 		Sgl_leftshift(opnd3,SGL_P-count);
239 		if (Sgl_iszero_hidden(opnd3)) {
240 			Sgl_leftshiftby1(opnd3);
241 			dest_exponent--;
242 		}
243 	}
244 	else {
245 		if (Sgl_iszero_hidden(opnd3)) {
246 			/* need to get one more bit of result */
247 			Sgl_leftshiftby1(opnd1);
248 			Sgl_leftshiftby1(opnd3);
249 			if (Sgl_iszero_sign(opnd1)) {
250 				Sgl_setone_lowmantissa(opnd3);
251 				Sgl_subtract(opnd1,opnd2,opnd1);
252 			}
253 			else Sgl_addition(opnd1,opnd2,opnd1);
254 			dest_exponent--;
255 		}
256 		if (Sgl_iszero_sign(opnd1)) guardbit = TRUE;
257 		stickybit = Sgl_all(opnd1);
258 	}
259 	inexact = guardbit | stickybit;
260 
261 	/*
262 	 * round result
263 	 */
264 	if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
265 		Sgl_clear_signexponent(opnd3);
266 		switch (Rounding_mode()) {
267 			case ROUNDPLUS:
268 				if (Sgl_iszero_sign(result))
269 					Sgl_increment_mantissa(opnd3);
270 				break;
271 			case ROUNDMINUS:
272 				if (Sgl_isone_sign(result))
273 					Sgl_increment_mantissa(opnd3);
274 				break;
275 			case ROUNDNEAREST:
276 				if (guardbit &&
277 				    (stickybit || Sgl_isone_lowmantissa(opnd3)))
278 					Sgl_increment_mantissa(opnd3);
279 		}
280 		if (Sgl_isone_hidden(opnd3)) dest_exponent++;
281 	}
282 	Sgl_set_mantissa(result,opnd3);
283 
284 	/*
285 	 * Test for overflow
286 	 */
287 	if (dest_exponent >= SGL_INFINITY_EXPONENT) {
288 		/* trap if OVERFLOWTRAP enabled */
289 		if (Is_overflowtrap_enabled()) {
290 			/*
291 			 * Adjust bias of result
292 			 */
293 			Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
294 			*dstptr = result;
295 			if (inexact) {
296 			    if (Is_inexacttrap_enabled())
297 				return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
298 			    else Set_inexactflag();
299 			}
300 			return(OVERFLOWEXCEPTION);
301 		}
302 		Set_overflowflag();
303 		/* set result to infinity or largest number */
304 		Sgl_setoverflow(result);
305 		inexact = TRUE;
306 	}
307 	/*
308 	 * Test for underflow
309 	 */
310 	else if (dest_exponent <= 0) {
311 		/* trap if UNDERFLOWTRAP enabled */
312 		if (Is_underflowtrap_enabled()) {
313 			/*
314 			 * Adjust bias of result
315 			 */
316 			Sgl_setwrapped_exponent(result,dest_exponent,unfl);
317 			*dstptr = result;
318 			if (inexact) {
319 			    if (Is_inexacttrap_enabled())
320 				return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
321 			    else Set_inexactflag();
322 			}
323 			return(UNDERFLOWEXCEPTION);
324 		}
325 
326 		/* Determine if should set underflow flag */
327 		is_tiny = TRUE;
328 		if (dest_exponent == 0 && inexact) {
329 			switch (Rounding_mode()) {
330 			case ROUNDPLUS:
331 				if (Sgl_iszero_sign(result)) {
332 					Sgl_increment(opnd3);
333 					if (Sgl_isone_hiddenoverflow(opnd3))
334 						is_tiny = FALSE;
335 					Sgl_decrement(opnd3);
336 				}
337 				break;
338 			case ROUNDMINUS:
339 				if (Sgl_isone_sign(result)) {
340 					Sgl_increment(opnd3);
341 					if (Sgl_isone_hiddenoverflow(opnd3))
342 						is_tiny = FALSE;
343 					Sgl_decrement(opnd3);
344 				}
345 				break;
346 			case ROUNDNEAREST:
347 				if (guardbit && (stickybit ||
348 				    Sgl_isone_lowmantissa(opnd3))) {
349 					Sgl_increment(opnd3);
350 					if (Sgl_isone_hiddenoverflow(opnd3))
351 						is_tiny = FALSE;
352 					Sgl_decrement(opnd3);
353 				}
354 				break;
355 			}
356 		}
357 
358 		/*
359 		 * denormalize result or set to signed zero
360 		 */
361 		stickybit = inexact;
362 		Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
363 
364 		/* return rounded number */
365 		if (inexact) {
366 			switch (Rounding_mode()) {
367 			case ROUNDPLUS:
368 				if (Sgl_iszero_sign(result)) {
369 					Sgl_increment(opnd3);
370 				}
371 				break;
372 			case ROUNDMINUS:
373 				if (Sgl_isone_sign(result))  {
374 					Sgl_increment(opnd3);
375 				}
376 				break;
377 			case ROUNDNEAREST:
378 				if (guardbit && (stickybit ||
379 				    Sgl_isone_lowmantissa(opnd3))) {
380 					Sgl_increment(opnd3);
381 				}
382 				break;
383 			}
384 			if (is_tiny) Set_underflowflag();
385 		}
386 		Sgl_set_exponentmantissa(result,opnd3);
387 	}
388 	else Sgl_set_exponent(result,dest_exponent);
389 	*dstptr = result;
390 	/* check for inexact */
391 	if (inexact) {
392 		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
393 		else  Set_inexactflag();
394 	}
395 	return(NOEXCEPTION);
396 }
397