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