1/*	$OpenBSD: dfrem.c,v 1.5 2002/05/07 22:19:30 mickey Exp $	*/
2/*
3  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
4  To anyone who acknowledges that this file is provided "AS IS"
5  without any express or implied warranty:
6      permission to use, copy, modify, and distribute this file
7  for any purpose is hereby granted without fee, provided that
8  the above copyright notice and this notice appears in all
9  copies, and that the name of Hewlett-Packard Company not be
10  used in advertising or publicity pertaining to distribution
11  of the software without specific, written prior permission.
12  Hewlett-Packard Company makes no representations about the
13  suitability of this software for any purpose.
14*/
15/* @(#)dfrem.c: Revision: 1.7.88.1 Date: 93/12/07 15:05:43 */
16
17#include "float.h"
18#include "dbl_float.h"
19
20/*
21 *  Double Precision Floating-point Remainder
22 */
23int
24dbl_frem(srcptr1,srcptr2,dstptr,status)
25	dbl_floating_point *srcptr1, *srcptr2, *dstptr;
26	unsigned int *status;
27{
28	register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
29	register unsigned int resultp1, resultp2;
30	register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
31	register int roundup = FALSE;
32
33	Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
34	Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
35	/*
36	 * check first operand for NaN's or infinity
37	 */
38	if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) {
39		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
40			if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
41				/* invalid since first operand is infinity */
42				if (Is_invalidtrap_enabled())
43					return(INVALIDEXCEPTION);
44				Set_invalidflag();
45				Dbl_makequietnan(resultp1,resultp2);
46				Dbl_copytoptr(resultp1,resultp2,dstptr);
47				return(NOEXCEPTION);
48			}
49		}
50		else {
51			/*
52			 * is NaN; signaling or quiet?
53			 */
54			if (Dbl_isone_signaling(opnd1p1)) {
55				/* trap if INVALIDTRAP enabled */
56				if (Is_invalidtrap_enabled())
57					return(INVALIDEXCEPTION);
58				/* make NaN quiet */
59				Set_invalidflag();
60				Dbl_set_quiet(opnd1p1);
61			}
62			/*
63			 * is second operand a signaling NaN?
64			 */
65			else if (Dbl_is_signalingnan(opnd2p1)) {
66				/* trap if INVALIDTRAP enabled */
67				if (Is_invalidtrap_enabled())
68					return(INVALIDEXCEPTION);
69				/* make NaN quiet */
70				Set_invalidflag();
71				Dbl_set_quiet(opnd2p1);
72				Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
73				return(NOEXCEPTION);
74			}
75			/*
76			 * return quiet NaN
77			 */
78			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
79			return(NOEXCEPTION);
80		}
81	}
82	/*
83	 * check second operand for NaN's or infinity
84	 */
85	if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
86		if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
87			/*
88			 * return first operand
89			 */
90			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
91			return(NOEXCEPTION);
92		}
93		/*
94		 * is NaN; signaling or quiet?
95		 */
96		if (Dbl_isone_signaling(opnd2p1)) {
97			/* trap if INVALIDTRAP enabled */
98			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
99			/* make NaN quiet */
100			Set_invalidflag();
101			Dbl_set_quiet(opnd2p1);
102		}
103		/*
104		 * return quiet NaN
105		 */
106		Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
107		return(NOEXCEPTION);
108	}
109	/*
110	 * check second operand for zero
111	 */
112	if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
113		/* invalid since second operand is zero */
114		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
115		Set_invalidflag();
116		Dbl_makequietnan(resultp1,resultp2);
117		Dbl_copytoptr(resultp1,resultp2,dstptr);
118		return(NOEXCEPTION);
119	}
120
121	/*
122	 * get sign of result
123	 */
124	resultp1 = opnd1p1;
125
126	/*
127	 * check for denormalized operands
128	 */
129	if (opnd1_exponent == 0) {
130		/* check for zero */
131		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
132			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
133			return(NOEXCEPTION);
134		}
135		/* normalize, then continue */
136		opnd1_exponent = 1;
137		Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
138	}
139	else {
140		Dbl_clear_signexponent_set_hidden(opnd1p1);
141	}
142	if (opnd2_exponent == 0) {
143		/* normalize, then continue */
144		opnd2_exponent = 1;
145		Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
146	}
147	else {
148		Dbl_clear_signexponent_set_hidden(opnd2p1);
149	}
150
151	/* find result exponent and divide step loop count */
152	dest_exponent = opnd2_exponent - 1;
153	stepcount = opnd1_exponent - opnd2_exponent;
154
155	/*
156	 * check for opnd1/opnd2 < 1
157	 */
158	if (stepcount < 0) {
159		/*
160		 * check for opnd1/opnd2 > 1/2
161		 *
162		 * In this case n will round to 1, so
163		 *    r = opnd1 - opnd2
164		 */
165		if (stepcount == -1 &&
166		    Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
167			/* set sign */
168			Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1);
169			/* align opnd2 with opnd1 */
170			Dbl_leftshiftby1(opnd2p1,opnd2p2);
171			Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,
172			 opnd2p1,opnd2p2);
173			/* now normalize */
174			while (Dbl_iszero_hidden(opnd2p1)) {
175				Dbl_leftshiftby1(opnd2p1,opnd2p2);
176				dest_exponent--;
177			}
178			Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
179			goto testforunderflow;
180		}
181		/*
182		 * opnd1/opnd2 <= 1/2
183		 *
184		 * In this case n will round to zero, so
185		 *    r = opnd1
186		 */
187		Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
188		dest_exponent = opnd1_exponent;
189		goto testforunderflow;
190	}
191
192	/*
193	 * Generate result
194	 *
195	 * Do iterative subtract until remainder is less than operand 2.
196	 */
197	while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) {
198		if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
199			Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
200		}
201		Dbl_leftshiftby1(opnd1p1,opnd1p2);
202	}
203	/*
204	 * Do last subtract, then determine which way to round if remainder
205	 * is exactly 1/2 of opnd2
206	 */
207	if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
208		Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
209		roundup = TRUE;
210	}
211	if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) {
212		/* division is exact, remainder is zero */
213		Dbl_setzero_exponentmantissa(resultp1,resultp2);
214		Dbl_copytoptr(resultp1,resultp2,dstptr);
215		return(NOEXCEPTION);
216	}
217
218	/*
219	 * Check for cases where opnd1/opnd2 < n
220	 *
221	 * In this case the result's sign will be opposite that of
222	 * opnd1.  The mantissa also needs some correction.
223	 */
224	Dbl_leftshiftby1(opnd1p1,opnd1p2);
225	if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
226		Dbl_invert_sign(resultp1);
227		Dbl_leftshiftby1(opnd2p1,opnd2p2);
228		Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2);
229	}
230	/* check for remainder being exactly 1/2 of opnd2 */
231	else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) {
232		Dbl_invert_sign(resultp1);
233	}
234
235	/* normalize result's mantissa */
236	while (Dbl_iszero_hidden(opnd1p1)) {
237		dest_exponent--;
238		Dbl_leftshiftby1(opnd1p1,opnd1p2);
239	}
240	Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
241
242	/*
243	 * Test for underflow
244	 */
245    testforunderflow:
246	if (dest_exponent <= 0) {
247		/* trap if UNDERFLOWTRAP enabled */
248		if (Is_underflowtrap_enabled()) {
249			/*
250			 * Adjust bias of result
251			 */
252			Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
253			/* frem is always exact */
254			Dbl_copytoptr(resultp1,resultp2,dstptr);
255			return(UNDERFLOWEXCEPTION);
256		}
257		/*
258		 * denormalize result or set to signed zero
259		 */
260		if (dest_exponent >= (1 - DBL_P)) {
261			Dbl_rightshift_exponentmantissa(resultp1,resultp2,
262			 1-dest_exponent);
263		}
264		else {
265			Dbl_setzero_exponentmantissa(resultp1,resultp2);
266		}
267	}
268	else Dbl_set_exponent(resultp1,dest_exponent);
269	Dbl_copytoptr(resultp1,resultp2,dstptr);
270	return(NOEXCEPTION);
271}
272