xref: /linux/arch/parisc/math-emu/dfdiv.c (revision ca64d84e93762f4e587e040a44ad9f6089afc777)
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/dfdiv.c		$Revision: 1.1 $
13  *
14  *  Purpose:
15  *	Double Precision Floating-point Divide
16  *
17  *  External Interfaces:
18  *	dbl_fdiv(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 "dbl_float.h"
31 
32 /*
33  *  Double Precision Floating-point Divide
34  */
35 
36 int
37 dbl_fdiv (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2,
38 	  dbl_floating_point * dstptr, unsigned int *status)
39 {
40 	register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
41 	register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
42 	register int dest_exponent, count;
43 	register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
44 	boolean is_tiny;
45 
46 	Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
47 	Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
48 	/*
49 	 * set sign bit of result
50 	 */
51 	if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
52 		Dbl_setnegativezerop1(resultp1);
53 	else Dbl_setzerop1(resultp1);
54 	/*
55 	 * check first operand for NaN's or infinity
56 	 */
57 	if (Dbl_isinfinity_exponent(opnd1p1)) {
58 		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
59 			if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
60 				if (Dbl_isinfinity(opnd2p1,opnd2p2)) {
61 					/*
62 					 * invalid since both operands
63 					 * are infinity
64 					 */
65 					if (Is_invalidtrap_enabled())
66                                 		return(INVALIDEXCEPTION);
67                                 	Set_invalidflag();
68                                 	Dbl_makequietnan(resultp1,resultp2);
69 					Dbl_copytoptr(resultp1,resultp2,dstptr);
70 					return(NOEXCEPTION);
71 				}
72 				/*
73 			 	 * return infinity
74 			 	 */
75 				Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
76 				Dbl_copytoptr(resultp1,resultp2,dstptr);
77 				return(NOEXCEPTION);
78 			}
79 		}
80 		else {
81                 	/*
82                  	 * is NaN; signaling or quiet?
83                  	 */
84                 	if (Dbl_isone_signaling(opnd1p1)) {
85                         	/* trap if INVALIDTRAP enabled */
86                         	if (Is_invalidtrap_enabled())
87                             		return(INVALIDEXCEPTION);
88                         	/* make NaN quiet */
89                         	Set_invalidflag();
90                         	Dbl_set_quiet(opnd1p1);
91                 	}
92 			/*
93 			 * is second operand a signaling NaN?
94 			 */
95 			else if (Dbl_is_signalingnan(opnd2p1)) {
96                         	/* trap if INVALIDTRAP enabled */
97                         	if (Is_invalidtrap_enabled())
98                             		return(INVALIDEXCEPTION);
99                         	/* make NaN quiet */
100                         	Set_invalidflag();
101                         	Dbl_set_quiet(opnd2p1);
102 				Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
103                 		return(NOEXCEPTION);
104 			}
105                 	/*
106                  	 * return quiet NaN
107                  	 */
108 			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
109                 	return(NOEXCEPTION);
110 		}
111 	}
112 	/*
113 	 * check second operand for NaN's or infinity
114 	 */
115 	if (Dbl_isinfinity_exponent(opnd2p1)) {
116 		if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
117 			/*
118 			 * return zero
119 			 */
120 			Dbl_setzero_exponentmantissa(resultp1,resultp2);
121 			Dbl_copytoptr(resultp1,resultp2,dstptr);
122 			return(NOEXCEPTION);
123 		}
124                 /*
125                  * is NaN; signaling or quiet?
126                  */
127                 if (Dbl_isone_signaling(opnd2p1)) {
128                         /* trap if INVALIDTRAP enabled */
129                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
130                         /* make NaN quiet */
131                         Set_invalidflag();
132                         Dbl_set_quiet(opnd2p1);
133                 }
134                 /*
135                  * return quiet NaN
136                  */
137 		Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
138                 return(NOEXCEPTION);
139 	}
140         /*
141          * check for division by zero
142          */
143         if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
144                 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
145                         /* invalid since both operands are zero */
146                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
147                         Set_invalidflag();
148                         Dbl_makequietnan(resultp1,resultp2);
149                         Dbl_copytoptr(resultp1,resultp2,dstptr);
150                         return(NOEXCEPTION);
151                 }
152                 if (Is_divisionbyzerotrap_enabled())
153                        	return(DIVISIONBYZEROEXCEPTION);
154                 Set_divisionbyzeroflag();
155                 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
156                 Dbl_copytoptr(resultp1,resultp2,dstptr);
157                 return(NOEXCEPTION);
158         }
159 	/*
160 	 * Generate exponent
161 	 */
162 	dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS;
163 
164 	/*
165 	 * Generate mantissa
166 	 */
167 	if (Dbl_isnotzero_exponent(opnd1p1)) {
168 		/* set hidden bit */
169 		Dbl_clear_signexponent_set_hidden(opnd1p1);
170 	}
171 	else {
172 		/* check for zero */
173 		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
174 			Dbl_setzero_exponentmantissa(resultp1,resultp2);
175 			Dbl_copytoptr(resultp1,resultp2,dstptr);
176 			return(NOEXCEPTION);
177 		}
178                 /* is denormalized, want to normalize */
179                 Dbl_clear_signexponent(opnd1p1);
180                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
181 		Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
182 	}
183 	/* opnd2 needs to have hidden bit set with msb in hidden bit */
184 	if (Dbl_isnotzero_exponent(opnd2p1)) {
185 		Dbl_clear_signexponent_set_hidden(opnd2p1);
186 	}
187 	else {
188                 /* is denormalized; want to normalize */
189                 Dbl_clear_signexponent(opnd2p1);
190                 Dbl_leftshiftby1(opnd2p1,opnd2p2);
191                 while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) {
192                         dest_exponent+=8;
193                         Dbl_leftshiftby8(opnd2p1,opnd2p2);
194                 }
195                 if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) {
196                         dest_exponent+=4;
197                         Dbl_leftshiftby4(opnd2p1,opnd2p2);
198                 }
199                 while (Dbl_iszero_hidden(opnd2p1)) {
200                         dest_exponent++;
201                         Dbl_leftshiftby1(opnd2p1,opnd2p2);
202                 }
203 	}
204 
205 	/* Divide the source mantissas */
206 
207 	/*
208 	 * A non-restoring divide algorithm is used.
209 	 */
210 	Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
211 	Dbl_setzero(opnd3p1,opnd3p2);
212 	for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) {
213 		Dbl_leftshiftby1(opnd1p1,opnd1p2);
214 		Dbl_leftshiftby1(opnd3p1,opnd3p2);
215 		if (Dbl_iszero_sign(opnd1p1)) {
216 			Dbl_setone_lowmantissap2(opnd3p2);
217 			Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
218 		}
219 		else {
220 			Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2);
221 		}
222 	}
223 	if (count <= DBL_P) {
224 		Dbl_leftshiftby1(opnd3p1,opnd3p2);
225 		Dbl_setone_lowmantissap2(opnd3p2);
226 		Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count));
227 		if (Dbl_iszero_hidden(opnd3p1)) {
228 			Dbl_leftshiftby1(opnd3p1,opnd3p2);
229 			dest_exponent--;
230 		}
231 	}
232 	else {
233 		if (Dbl_iszero_hidden(opnd3p1)) {
234 			/* need to get one more bit of result */
235 			Dbl_leftshiftby1(opnd1p1,opnd1p2);
236 			Dbl_leftshiftby1(opnd3p1,opnd3p2);
237 			if (Dbl_iszero_sign(opnd1p1)) {
238 				Dbl_setone_lowmantissap2(opnd3p2);
239 				Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
240 			}
241 			else {
242 				Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
243 			}
244 			dest_exponent--;
245 		}
246 		if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE;
247 		stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2);
248 	}
249 	inexact = guardbit | stickybit;
250 
251 	/*
252 	 * round result
253 	 */
254 	if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
255 		Dbl_clear_signexponent(opnd3p1);
256 		switch (Rounding_mode()) {
257 			case ROUNDPLUS:
258 				if (Dbl_iszero_sign(resultp1))
259 					Dbl_increment(opnd3p1,opnd3p2);
260 				break;
261 			case ROUNDMINUS:
262 				if (Dbl_isone_sign(resultp1))
263 					Dbl_increment(opnd3p1,opnd3p2);
264 				break;
265 			case ROUNDNEAREST:
266 				if (guardbit && (stickybit ||
267 				    Dbl_isone_lowmantissap2(opnd3p2))) {
268 			      		Dbl_increment(opnd3p1,opnd3p2);
269 				}
270 		}
271 		if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
272 	}
273 	Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);
274 
275         /*
276          * Test for overflow
277          */
278 	if (dest_exponent >= DBL_INFINITY_EXPONENT) {
279                 /* trap if OVERFLOWTRAP enabled */
280                 if (Is_overflowtrap_enabled()) {
281                         /*
282                          * Adjust bias of result
283                          */
284                         Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
285                         Dbl_copytoptr(resultp1,resultp2,dstptr);
286                         if (inexact)
287                             if (Is_inexacttrap_enabled())
288                                 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
289                             else Set_inexactflag();
290                         return(OVERFLOWEXCEPTION);
291                 }
292 		Set_overflowflag();
293                 /* set result to infinity or largest number */
294 		Dbl_setoverflow(resultp1,resultp2);
295 		inexact = TRUE;
296 	}
297         /*
298          * Test for underflow
299          */
300 	else if (dest_exponent <= 0) {
301                 /* trap if UNDERFLOWTRAP enabled */
302                 if (Is_underflowtrap_enabled()) {
303                         /*
304                          * Adjust bias of result
305                          */
306                         Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
307                         Dbl_copytoptr(resultp1,resultp2,dstptr);
308                         if (inexact)
309                             if (Is_inexacttrap_enabled())
310                                 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
311                             else Set_inexactflag();
312                         return(UNDERFLOWEXCEPTION);
313                 }
314 
315 		/* Determine if should set underflow flag */
316 		is_tiny = TRUE;
317 		if (dest_exponent == 0 && inexact) {
318 			switch (Rounding_mode()) {
319 			case ROUNDPLUS:
320 				if (Dbl_iszero_sign(resultp1)) {
321 					Dbl_increment(opnd3p1,opnd3p2);
322 					if (Dbl_isone_hiddenoverflow(opnd3p1))
323                 			    is_tiny = FALSE;
324 					Dbl_decrement(opnd3p1,opnd3p2);
325 				}
326 				break;
327 			case ROUNDMINUS:
328 				if (Dbl_isone_sign(resultp1)) {
329 					Dbl_increment(opnd3p1,opnd3p2);
330 					if (Dbl_isone_hiddenoverflow(opnd3p1))
331                 			    is_tiny = FALSE;
332 					Dbl_decrement(opnd3p1,opnd3p2);
333 				}
334 				break;
335 			case ROUNDNEAREST:
336 				if (guardbit && (stickybit ||
337 				    Dbl_isone_lowmantissap2(opnd3p2))) {
338 				      	Dbl_increment(opnd3p1,opnd3p2);
339 					if (Dbl_isone_hiddenoverflow(opnd3p1))
340                 			    is_tiny = FALSE;
341 					Dbl_decrement(opnd3p1,opnd3p2);
342 				}
343 				break;
344 			}
345 		}
346 
347                 /*
348                  * denormalize result or set to signed zero
349                  */
350 		stickybit = inexact;
351 		Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
352 		 stickybit,inexact);
353 
354 		/* return rounded number */
355 		if (inexact) {
356 			switch (Rounding_mode()) {
357 			case ROUNDPLUS:
358 				if (Dbl_iszero_sign(resultp1)) {
359 					Dbl_increment(opnd3p1,opnd3p2);
360 				}
361 				break;
362 			case ROUNDMINUS:
363 				if (Dbl_isone_sign(resultp1)) {
364 					Dbl_increment(opnd3p1,opnd3p2);
365 				}
366 				break;
367 			case ROUNDNEAREST:
368 				if (guardbit && (stickybit ||
369 				    Dbl_isone_lowmantissap2(opnd3p2))) {
370 			      		Dbl_increment(opnd3p1,opnd3p2);
371 				}
372 				break;
373 			}
374                 	if (is_tiny) Set_underflowflag();
375                 }
376 		Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
377 	}
378 	else Dbl_set_exponent(resultp1,dest_exponent);
379 	Dbl_copytoptr(resultp1,resultp2,dstptr);
380 
381 	/* check for inexact */
382 	if (inexact) {
383 		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
384 		else Set_inexactflag();
385 	}
386 	return(NOEXCEPTION);
387 }
388