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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
20 */
21/*
22 * BEGIN_DESC
23 *
24 *  File:
25 *	@(#)	pa/spmath/dfrem.c		$Revision: 1.1.1.1 $
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
50int
51dbl_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