1 /* $NetBSD: dfrem.c,v 1.5 2012/02/04 17:03:09 skrll Exp $ */
2
3 /* $OpenBSD: dfrem.c,v 1.4 2001/03/29 03:58:17 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: dfrem.c,v 1.5 2012/02/04 17:03:09 skrll Exp $");
46
47 #include "../spmath/float.h"
48 #include "../spmath/dbl_float.h"
49
50 /*
51 * Double Precision Floating-point Remainder
52 */
53 int
dbl_frem(dbl_floating_point * srcptr1,dbl_floating_point * srcptr2,dbl_floating_point * dstptr,unsigned int * status)54 dbl_frem(dbl_floating_point *srcptr1, dbl_floating_point *srcptr2,
55 dbl_floating_point *dstptr, unsigned int *status)
56 {
57 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
58 register unsigned int resultp1, resultp2;
59 register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
60 register int roundup = false;
61
62 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
63 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
64 /*
65 * check first operand for NaN's or infinity
66 */
67 if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) {
68 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
69 if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
70 /* invalid since first operand is infinity */
71 if (Is_invalidtrap_enabled())
72 return(INVALIDEXCEPTION);
73 Set_invalidflag();
74 Dbl_makequietnan(resultp1,resultp2);
75 Dbl_copytoptr(resultp1,resultp2,dstptr);
76 return(NOEXCEPTION);
77 }
78 }
79 else {
80 /*
81 * is NaN; signaling or quiet?
82 */
83 if (Dbl_isone_signaling(opnd1p1)) {
84 /* trap if INVALIDTRAP enabled */
85 if (Is_invalidtrap_enabled())
86 return(INVALIDEXCEPTION);
87 /* make NaN quiet */
88 Set_invalidflag();
89 Dbl_set_quiet(opnd1p1);
90 }
91 /*
92 * is second operand a signaling NaN?
93 */
94 else if (Dbl_is_signalingnan(opnd2p1)) {
95 /* trap if INVALIDTRAP enabled */
96 if (Is_invalidtrap_enabled())
97 return(INVALIDEXCEPTION);
98 /* make NaN quiet */
99 Set_invalidflag();
100 Dbl_set_quiet(opnd2p1);
101 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
102 return(NOEXCEPTION);
103 }
104 /*
105 * return quiet NaN
106 */
107 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
108 return(NOEXCEPTION);
109 }
110 }
111 /*
112 * check second operand for NaN's or infinity
113 */
114 if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
115 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
116 /*
117 * return first operand
118 */
119 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
120 return(NOEXCEPTION);
121 }
122 /*
123 * is NaN; signaling or quiet?
124 */
125 if (Dbl_isone_signaling(opnd2p1)) {
126 /* trap if INVALIDTRAP enabled */
127 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
128 /* make NaN quiet */
129 Set_invalidflag();
130 Dbl_set_quiet(opnd2p1);
131 }
132 /*
133 * return quiet NaN
134 */
135 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
136 return(NOEXCEPTION);
137 }
138 /*
139 * check second operand for zero
140 */
141 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
142 /* invalid since second operand is zero */
143 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
144 Set_invalidflag();
145 Dbl_makequietnan(resultp1,resultp2);
146 Dbl_copytoptr(resultp1,resultp2,dstptr);
147 return(NOEXCEPTION);
148 }
149
150 /*
151 * get sign of result
152 */
153 resultp1 = opnd1p1;
154
155 /*
156 * check for denormalized operands
157 */
158 if (opnd1_exponent == 0) {
159 /* check for zero */
160 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
161 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
162 return(NOEXCEPTION);
163 }
164 /* normalize, then continue */
165 opnd1_exponent = 1;
166 Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
167 }
168 else {
169 Dbl_clear_signexponent_set_hidden(opnd1p1);
170 }
171 if (opnd2_exponent == 0) {
172 /* normalize, then continue */
173 opnd2_exponent = 1;
174 Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
175 }
176 else {
177 Dbl_clear_signexponent_set_hidden(opnd2p1);
178 }
179
180 /* find result exponent and divide step loop count */
181 dest_exponent = opnd2_exponent - 1;
182 stepcount = opnd1_exponent - opnd2_exponent;
183
184 /*
185 * check for opnd1/opnd2 < 1
186 */
187 if (stepcount < 0) {
188 /*
189 * check for opnd1/opnd2 > 1/2
190 *
191 * In this case n will round to 1, so
192 * r = opnd1 - opnd2
193 */
194 if (stepcount == -1 &&
195 Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
196 /* set sign */
197 Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1);
198 /* align opnd2 with opnd1 */
199 Dbl_leftshiftby1(opnd2p1,opnd2p2);
200 Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,
201 opnd2p1,opnd2p2);
202 /* now normalize */
203 while (Dbl_iszero_hidden(opnd2p1)) {
204 Dbl_leftshiftby1(opnd2p1,opnd2p2);
205 dest_exponent--;
206 }
207 Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
208 goto testforunderflow;
209 }
210 /*
211 * opnd1/opnd2 <= 1/2
212 *
213 * In this case n will round to zero, so
214 * r = opnd1
215 */
216 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
217 dest_exponent = opnd1_exponent;
218 goto testforunderflow;
219 }
220
221 /*
222 * Generate result
223 *
224 * Do iterative subtract until remainder is less than operand 2.
225 */
226 while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) {
227 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
228 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
229 }
230 Dbl_leftshiftby1(opnd1p1,opnd1p2);
231 }
232 /*
233 * Do last subtract, then determine which way to round if remainder
234 * is exactly 1/2 of opnd2
235 */
236 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
237 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
238 roundup = true;
239 }
240 if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) {
241 /* division is exact, remainder is zero */
242 Dbl_setzero_exponentmantissa(resultp1,resultp2);
243 Dbl_copytoptr(resultp1,resultp2,dstptr);
244 return(NOEXCEPTION);
245 }
246
247 /*
248 * Check for cases where opnd1/opnd2 < n
249 *
250 * In this case the result's sign will be opposite that of
251 * opnd1. The mantissa also needs some correction.
252 */
253 Dbl_leftshiftby1(opnd1p1,opnd1p2);
254 if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
255 Dbl_invert_sign(resultp1);
256 Dbl_leftshiftby1(opnd2p1,opnd2p2);
257 Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2);
258 }
259 /* check for remainder being exactly 1/2 of opnd2 */
260 else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) {
261 Dbl_invert_sign(resultp1);
262 }
263
264 /* normalize result's mantissa */
265 while (Dbl_iszero_hidden(opnd1p1)) {
266 dest_exponent--;
267 Dbl_leftshiftby1(opnd1p1,opnd1p2);
268 }
269 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
270
271 /*
272 * Test for underflow
273 */
274 testforunderflow:
275 if (dest_exponent <= 0) {
276 /* trap if UNDERFLOWTRAP enabled */
277 if (Is_underflowtrap_enabled()) {
278 /*
279 * Adjust bias of result
280 */
281 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
282 /* frem is always exact */
283 Dbl_copytoptr(resultp1,resultp2,dstptr);
284 return(UNDERFLOWEXCEPTION);
285 }
286 /*
287 * denormalize result or set to signed zero
288 */
289 if (dest_exponent >= (1 - DBL_P)) {
290 Dbl_rightshift_exponentmantissa(resultp1,resultp2,
291 1-dest_exponent);
292 }
293 else {
294 Dbl_setzero_exponentmantissa(resultp1,resultp2);
295 }
296 }
297 else Dbl_set_exponent(resultp1,dest_exponent);
298 Dbl_copytoptr(resultp1,resultp2,dstptr);
299 return(NOEXCEPTION);
300 }
301