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 */
23 int
dbl_frem(srcptr1,srcptr2,dstptr,status)24 dbl_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