xref: /netbsd-src/sys/arch/hppa/spmath/sfdiv.c (revision b1c86f5f087524e68db12794ee9c3e3da1ab17a0)
1 /*	$NetBSD: sfdiv.c,v 1.4 2007/02/22 05:46:30 thorpej 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 <sys/cdefs.h>
45 __KERNEL_RCSID(0, "$NetBSD: sfdiv.c,v 1.4 2007/02/22 05:46:30 thorpej Exp $");
46 
47 #include "../spmath/float.h"
48 #include "../spmath/sgl_float.h"
49 
50 /*
51  *  Single Precision Floating-point Divide
52  */
53 int
54 sgl_fdiv(srcptr1,srcptr2,dstptr,status)
55 
56 sgl_floating_point *srcptr1, *srcptr2, *dstptr;
57 unsigned int *status;
58 {
59 	register unsigned int opnd1, opnd2, opnd3, result;
60 	register int dest_exponent, count;
61 	register int inexact = false, guardbit = false, stickybit = false;
62 	int is_tiny;
63 
64 	opnd1 = *srcptr1;
65 	opnd2 = *srcptr2;
66 	/*
67 	 * set sign bit of result
68 	 */
69 	if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);
70 	else Sgl_setzero(result);
71 	/*
72 	 * check first operand for NaN's or infinity
73 	 */
74 	if (Sgl_isinfinity_exponent(opnd1)) {
75 		if (Sgl_iszero_mantissa(opnd1)) {
76 			if (Sgl_isnotnan(opnd2)) {
77 				if (Sgl_isinfinity(opnd2)) {
78 					/*
79 					 * invalid since both operands
80 					 * are infinity
81 					 */
82 					if (Is_invalidtrap_enabled())
83 						return(INVALIDEXCEPTION);
84 					Set_invalidflag();
85 					Sgl_makequietnan(result);
86 					*dstptr = result;
87 					return(NOEXCEPTION);
88 				}
89 				/*
90 				 * return infinity
91 				 */
92 				Sgl_setinfinity_exponentmantissa(result);
93 				*dstptr = result;
94 				return(NOEXCEPTION);
95 			}
96 		}
97 		else {
98 			/*
99 			 * is NaN; signaling or quiet?
100 			 */
101 			if (Sgl_isone_signaling(opnd1)) {
102 				/* trap if INVALIDTRAP enabled */
103 				if (Is_invalidtrap_enabled())
104 					return(INVALIDEXCEPTION);
105 				/* make NaN quiet */
106 				Set_invalidflag();
107 				Sgl_set_quiet(opnd1);
108 			}
109 			/*
110 			 * is second operand a signaling NaN?
111 			 */
112 			else if (Sgl_is_signalingnan(opnd2)) {
113 				/* trap if INVALIDTRAP enabled */
114 				if (Is_invalidtrap_enabled())
115 					return(INVALIDEXCEPTION);
116 				/* make NaN quiet */
117 				Set_invalidflag();
118 				Sgl_set_quiet(opnd2);
119 				*dstptr = opnd2;
120 				return(NOEXCEPTION);
121 			}
122 			/*
123 			 * return quiet NaN
124 			 */
125 			*dstptr = opnd1;
126 			return(NOEXCEPTION);
127 		}
128 	}
129 	/*
130 	 * check second operand for NaN's or infinity
131 	 */
132 	if (Sgl_isinfinity_exponent(opnd2)) {
133 		if (Sgl_iszero_mantissa(opnd2)) {
134 			/*
135 			 * return zero
136 			 */
137 			Sgl_setzero_exponentmantissa(result);
138 			*dstptr = result;
139 			return(NOEXCEPTION);
140 		}
141 		/*
142 		 * is NaN; signaling or quiet?
143 		 */
144 		if (Sgl_isone_signaling(opnd2)) {
145 			/* trap if INVALIDTRAP enabled */
146 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
147 			/* make NaN quiet */
148 			Set_invalidflag();
149 			Sgl_set_quiet(opnd2);
150 		}
151 		/*
152 		 * return quiet NaN
153 		 */
154 		*dstptr = opnd2;
155 		return(NOEXCEPTION);
156 	}
157 	/*
158 	 * check for division by zero
159 	 */
160 	if (Sgl_iszero_exponentmantissa(opnd2)) {
161 		if (Sgl_iszero_exponentmantissa(opnd1)) {
162 			/* invalid since both operands are zero */
163 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
164 			Set_invalidflag();
165 			Sgl_makequietnan(result);
166 			*dstptr = result;
167 			return(NOEXCEPTION);
168 		}
169 		if (Is_divisionbyzerotrap_enabled())
170 			return(DIVISIONBYZEROEXCEPTION);
171 		Set_divisionbyzeroflag();
172 		Sgl_setinfinity_exponentmantissa(result);
173 		*dstptr = result;
174 		return(NOEXCEPTION);
175 	}
176 	/*
177 	 * Generate exponent
178 	 */
179 	dest_exponent = Sgl_exponent(opnd1) - Sgl_exponent(opnd2) + SGL_BIAS;
180 
181 	/*
182 	 * Generate mantissa
183 	 */
184 	if (Sgl_isnotzero_exponent(opnd1)) {
185 		/* set hidden bit */
186 		Sgl_clear_signexponent_set_hidden(opnd1);
187 	}
188 	else {
189 		/* check for zero */
190 		if (Sgl_iszero_mantissa(opnd1)) {
191 			Sgl_setzero_exponentmantissa(result);
192 			*dstptr = result;
193 			return(NOEXCEPTION);
194 		}
195 		/* is denormalized; want to normalize */
196 		Sgl_clear_signexponent(opnd1);
197 		Sgl_leftshiftby1(opnd1);
198 		Sgl_normalize(opnd1,dest_exponent);
199 	}
200 	/* opnd2 needs to have hidden bit set with msb in hidden bit */
201 	if (Sgl_isnotzero_exponent(opnd2)) {
202 		Sgl_clear_signexponent_set_hidden(opnd2);
203 	}
204 	else {
205 		/* is denormalized; want to normalize */
206 		Sgl_clear_signexponent(opnd2);
207 		Sgl_leftshiftby1(opnd2);
208 		while(Sgl_iszero_hiddenhigh7mantissa(opnd2)) {
209 			Sgl_leftshiftby8(opnd2);
210 			dest_exponent += 8;
211 		}
212 		if(Sgl_iszero_hiddenhigh3mantissa(opnd2)) {
213 			Sgl_leftshiftby4(opnd2);
214 			dest_exponent += 4;
215 		}
216 		while(Sgl_iszero_hidden(opnd2)) {
217 			Sgl_leftshiftby1(opnd2);
218 			dest_exponent += 1;
219 		}
220 	}
221 
222 	/* Divide the source mantissas */
223 
224 	/*
225 	 * A non_restoring divide algorithm is used.
226 	 */
227 	Sgl_subtract(opnd1,opnd2,opnd1);
228 	Sgl_setzero(opnd3);
229 	for (count=1;count<=SGL_P && Sgl_all(opnd1);count++) {
230 		Sgl_leftshiftby1(opnd1);
231 		Sgl_leftshiftby1(opnd3);
232 		if (Sgl_iszero_sign(opnd1)) {
233 			Sgl_setone_lowmantissa(opnd3);
234 			Sgl_subtract(opnd1,opnd2,opnd1);
235 		}
236 		else Sgl_addition(opnd1,opnd2,opnd1);
237 	}
238 	if (count <= SGL_P) {
239 		Sgl_leftshiftby1(opnd3);
240 		Sgl_setone_lowmantissa(opnd3);
241 		Sgl_leftshift(opnd3,SGL_P-count);
242 		if (Sgl_iszero_hidden(opnd3)) {
243 			Sgl_leftshiftby1(opnd3);
244 			dest_exponent--;
245 		}
246 	}
247 	else {
248 		if (Sgl_iszero_hidden(opnd3)) {
249 			/* need to get one more bit of result */
250 			Sgl_leftshiftby1(opnd1);
251 			Sgl_leftshiftby1(opnd3);
252 			if (Sgl_iszero_sign(opnd1)) {
253 				Sgl_setone_lowmantissa(opnd3);
254 				Sgl_subtract(opnd1,opnd2,opnd1);
255 			}
256 			else Sgl_addition(opnd1,opnd2,opnd1);
257 			dest_exponent--;
258 		}
259 		if (Sgl_iszero_sign(opnd1)) guardbit = true;
260 		stickybit = Sgl_all(opnd1);
261 	}
262 	inexact = guardbit | stickybit;
263 
264 	/*
265 	 * round result
266 	 */
267 	if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
268 		Sgl_clear_signexponent(opnd3);
269 		switch (Rounding_mode()) {
270 			case ROUNDPLUS:
271 				if (Sgl_iszero_sign(result))
272 					Sgl_increment_mantissa(opnd3);
273 				break;
274 			case ROUNDMINUS:
275 				if (Sgl_isone_sign(result))
276 					Sgl_increment_mantissa(opnd3);
277 				break;
278 			case ROUNDNEAREST:
279 				if (guardbit &&
280 				    (stickybit || Sgl_isone_lowmantissa(opnd3)))
281 					Sgl_increment_mantissa(opnd3);
282 		}
283 		if (Sgl_isone_hidden(opnd3)) dest_exponent++;
284 	}
285 	Sgl_set_mantissa(result,opnd3);
286 
287 	/*
288 	 * Test for overflow
289 	 */
290 	if (dest_exponent >= SGL_INFINITY_EXPONENT) {
291 		/* trap if OVERFLOWTRAP enabled */
292 		if (Is_overflowtrap_enabled()) {
293 			/*
294 			 * Adjust bias of result
295 			 */
296 			Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
297 			*dstptr = result;
298 			if (inexact) {
299 			    if (Is_inexacttrap_enabled())
300 				return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
301 			    else Set_inexactflag();
302 			}
303 			return(OVERFLOWEXCEPTION);
304 		}
305 		Set_overflowflag();
306 		/* set result to infinity or largest number */
307 		Sgl_setoverflow(result);
308 		inexact = true;
309 	}
310 	/*
311 	 * Test for underflow
312 	 */
313 	else if (dest_exponent <= 0) {
314 		/* trap if UNDERFLOWTRAP enabled */
315 		if (Is_underflowtrap_enabled()) {
316 			/*
317 			 * Adjust bias of result
318 			 */
319 			Sgl_setwrapped_exponent(result,dest_exponent,unfl);
320 			*dstptr = result;
321 			if (inexact) {
322 			    if (Is_inexacttrap_enabled())
323 				return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
324 			    else 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 (Sgl_iszero_sign(result)) {
335 					Sgl_increment(opnd3);
336 					if (Sgl_isone_hiddenoverflow(opnd3))
337 						is_tiny = false;
338 					Sgl_decrement(opnd3);
339 				}
340 				break;
341 			case ROUNDMINUS:
342 				if (Sgl_isone_sign(result)) {
343 					Sgl_increment(opnd3);
344 					if (Sgl_isone_hiddenoverflow(opnd3))
345 						is_tiny = false;
346 					Sgl_decrement(opnd3);
347 				}
348 				break;
349 			case ROUNDNEAREST:
350 				if (guardbit && (stickybit ||
351 				    Sgl_isone_lowmantissa(opnd3))) {
352 					Sgl_increment(opnd3);
353 					if (Sgl_isone_hiddenoverflow(opnd3))
354 						is_tiny = false;
355 					Sgl_decrement(opnd3);
356 				}
357 				break;
358 			}
359 		}
360 
361 		/*
362 		 * denormalize result or set to signed zero
363 		 */
364 		stickybit = inexact;
365 		Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
366 
367 		/* return rounded number */
368 		if (inexact) {
369 			switch (Rounding_mode()) {
370 			case ROUNDPLUS:
371 				if (Sgl_iszero_sign(result)) {
372 					Sgl_increment(opnd3);
373 				}
374 				break;
375 			case ROUNDMINUS:
376 				if (Sgl_isone_sign(result))  {
377 					Sgl_increment(opnd3);
378 				}
379 				break;
380 			case ROUNDNEAREST:
381 				if (guardbit && (stickybit ||
382 				    Sgl_isone_lowmantissa(opnd3))) {
383 					Sgl_increment(opnd3);
384 				}
385 				break;
386 			}
387 			if (is_tiny) Set_underflowflag();
388 		}
389 		Sgl_set_exponentmantissa(result,opnd3);
390 	}
391 	else Sgl_set_exponent(result,dest_exponent);
392 	*dstptr = result;
393 	/* check for inexact */
394 	if (inexact) {
395 		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
396 		else  Set_inexactflag();
397 	}
398 	return(NOEXCEPTION);
399 }
400