1 /*
2  * Linux/PA-RISC Project (http://www.parisc-linux.org/)
3  *
4  * Floating-point emulation code
5  *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
6  *
7  *    This program is free software; you can redistribute it and/or modify
8  *    it under the terms of the GNU General Public License as published by
9  *    the Free Software Foundation; either version 2, or (at your option)
10  *    any later version.
11  *
12  *    This program is distributed in the hope that it will be useful,
13  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  *    GNU General Public License for more details.
16  *
17  *    You should have received a copy of the GNU General Public License
18  *    along with this program; if not, write to the Free Software
19  *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
20  */
21 /*
22  * BEGIN_DESC
23  *
24  *  File:
25  *	@(#)	pa/spmath/dfrem.c		$Revision: $
26  *
27  *  Purpose:
28  *	Double Precision Floating-point Remainder
29  *
30  *  External Interfaces:
31  *	dbl_frem(srcptr1,srcptr2,dstptr,status)
32  *
33  *  Internal Interfaces:
34  *
35  *  Theory:
36  *	<<please update with a overview of the operation of this file>>
37  *
38  * END_DESC
39 */
40 
41 
42 
43 #include "float.h"
44 #include "dbl_float.h"
45 
46 /*
47  *  Double Precision Floating-point Remainder
48  */
49 
50 int
dbl_frem(dbl_floating_point * srcptr1,dbl_floating_point * srcptr2,dbl_floating_point * dstptr,unsigned int * status)51 dbl_frem (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2,
52 	  dbl_floating_point * dstptr, unsigned int *status)
53 {
54 	register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
55 	register unsigned int resultp1, resultp2;
56 	register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
57 	register boolean roundup = FALSE;
58 
59 	Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
60 	Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
61 	/*
62 	 * check first operand for NaN's or infinity
63 	 */
64 	if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) {
65 		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
66 			if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
67 				/* invalid since first operand is infinity */
68 				if (Is_invalidtrap_enabled())
69                                 	return(INVALIDEXCEPTION);
70                                 Set_invalidflag();
71                                 Dbl_makequietnan(resultp1,resultp2);
72 				Dbl_copytoptr(resultp1,resultp2,dstptr);
73 				return(NOEXCEPTION);
74 			}
75 		}
76 		else {
77                 	/*
78                  	 * is NaN; signaling or quiet?
79                  	 */
80                 	if (Dbl_isone_signaling(opnd1p1)) {
81                         	/* trap if INVALIDTRAP enabled */
82                         	if (Is_invalidtrap_enabled())
83                             		return(INVALIDEXCEPTION);
84                         	/* make NaN quiet */
85                         	Set_invalidflag();
86                         	Dbl_set_quiet(opnd1p1);
87                 	}
88 			/*
89 			 * is second operand a signaling NaN?
90 			 */
91 			else if (Dbl_is_signalingnan(opnd2p1)) {
92                         	/* trap if INVALIDTRAP enabled */
93                         	if (Is_invalidtrap_enabled())
94                             		return(INVALIDEXCEPTION);
95                         	/* make NaN quiet */
96                         	Set_invalidflag();
97                         	Dbl_set_quiet(opnd2p1);
98 				Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
99                 		return(NOEXCEPTION);
100 			}
101                 	/*
102                  	 * return quiet NaN
103                  	 */
104 			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
105                 	return(NOEXCEPTION);
106 		}
107 	}
108 	/*
109 	 * check second operand for NaN's or infinity
110 	 */
111 	if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
112 		if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
113 			/*
114 			 * return first operand
115 			 */
116 			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
117 			return(NOEXCEPTION);
118 		}
119                 /*
120                  * is NaN; signaling or quiet?
121                  */
122                 if (Dbl_isone_signaling(opnd2p1)) {
123                         /* trap if INVALIDTRAP enabled */
124                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
125                         /* make NaN quiet */
126                         Set_invalidflag();
127                         Dbl_set_quiet(opnd2p1);
128                 }
129                 /*
130                  * return quiet NaN
131                  */
132 		Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
133                 return(NOEXCEPTION);
134 	}
135 	/*
136 	 * check second operand for zero
137 	 */
138 	if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
139 		/* invalid since second operand is zero */
140 		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
141                 Set_invalidflag();
142                 Dbl_makequietnan(resultp1,resultp2);
143 		Dbl_copytoptr(resultp1,resultp2,dstptr);
144 		return(NOEXCEPTION);
145 	}
146 
147 	/*
148 	 * get sign of result
149 	 */
150 	resultp1 = opnd1p1;
151 
152 	/*
153 	 * check for denormalized operands
154 	 */
155 	if (opnd1_exponent == 0) {
156 		/* check for zero */
157 		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
158 			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
159 			return(NOEXCEPTION);
160 		}
161 		/* normalize, then continue */
162 		opnd1_exponent = 1;
163 		Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
164 	}
165 	else {
166 		Dbl_clear_signexponent_set_hidden(opnd1p1);
167 	}
168 	if (opnd2_exponent == 0) {
169 		/* normalize, then continue */
170 		opnd2_exponent = 1;
171 		Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
172 	}
173 	else {
174 		Dbl_clear_signexponent_set_hidden(opnd2p1);
175 	}
176 
177 	/* find result exponent and divide step loop count */
178 	dest_exponent = opnd2_exponent - 1;
179 	stepcount = opnd1_exponent - opnd2_exponent;
180 
181 	/*
182 	 * check for opnd1/opnd2 < 1
183 	 */
184 	if (stepcount < 0) {
185 		/*
186 		 * check for opnd1/opnd2 > 1/2
187 		 *
188 		 * In this case n will round to 1, so
189 		 *    r = opnd1 - opnd2
190 		 */
191 		if (stepcount == -1 &&
192 		    Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
193 			/* set sign */
194 			Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1);
195 			/* align opnd2 with opnd1 */
196 			Dbl_leftshiftby1(opnd2p1,opnd2p2);
197 			Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,
198 			 opnd2p1,opnd2p2);
199 			/* now normalize */
200                 	while (Dbl_iszero_hidden(opnd2p1)) {
201                         	Dbl_leftshiftby1(opnd2p1,opnd2p2);
202                         	dest_exponent--;
203 			}
204 			Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
205 			goto testforunderflow;
206 		}
207 		/*
208 		 * opnd1/opnd2 <= 1/2
209 		 *
210 		 * In this case n will round to zero, so
211 		 *    r = opnd1
212 		 */
213 		Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
214 		dest_exponent = opnd1_exponent;
215 		goto testforunderflow;
216 	}
217 
218 	/*
219 	 * Generate result
220 	 *
221 	 * Do iterative subtract until remainder is less than operand 2.
222 	 */
223 	while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) {
224 		if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
225 			Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
226 		}
227 		Dbl_leftshiftby1(opnd1p1,opnd1p2);
228 	}
229 	/*
230 	 * Do last subtract, then determine which way to round if remainder
231 	 * is exactly 1/2 of opnd2
232 	 */
233 	if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
234 		Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
235 		roundup = TRUE;
236 	}
237 	if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) {
238 		/* division is exact, remainder is zero */
239 		Dbl_setzero_exponentmantissa(resultp1,resultp2);
240 		Dbl_copytoptr(resultp1,resultp2,dstptr);
241 		return(NOEXCEPTION);
242 	}
243 
244 	/*
245 	 * Check for cases where opnd1/opnd2 < n
246 	 *
247 	 * In this case the result's sign will be opposite that of
248 	 * opnd1.  The mantissa also needs some correction.
249 	 */
250 	Dbl_leftshiftby1(opnd1p1,opnd1p2);
251 	if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
252 		Dbl_invert_sign(resultp1);
253 		Dbl_leftshiftby1(opnd2p1,opnd2p2);
254 		Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2);
255 	}
256 	/* check for remainder being exactly 1/2 of opnd2 */
257 	else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) {
258 		Dbl_invert_sign(resultp1);
259 	}
260 
261 	/* normalize result's mantissa */
262         while (Dbl_iszero_hidden(opnd1p1)) {
263                 dest_exponent--;
264                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
265         }
266 	Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
267 
268         /*
269          * Test for underflow
270          */
271     testforunderflow:
272 	if (dest_exponent <= 0) {
273                 /* trap if UNDERFLOWTRAP enabled */
274                 if (Is_underflowtrap_enabled()) {
275                         /*
276                          * Adjust bias of result
277                          */
278                         Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
279 			/* frem is always exact */
280 			Dbl_copytoptr(resultp1,resultp2,dstptr);
281 			return(UNDERFLOWEXCEPTION);
282                 }
283                 /*
284                  * denormalize result or set to signed zero
285                  */
286                 if (dest_exponent >= (1 - DBL_P)) {
287 			Dbl_rightshift_exponentmantissa(resultp1,resultp2,
288 			 1-dest_exponent);
289                 }
290                 else {
291 			Dbl_setzero_exponentmantissa(resultp1,resultp2);
292 		}
293 	}
294 	else Dbl_set_exponent(resultp1,dest_exponent);
295 	Dbl_copytoptr(resultp1,resultp2,dstptr);
296 	return(NOEXCEPTION);
297 }
298