Commit | Line | Data |
---|---|---|
1da177e4 LT |
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 $ | |
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 | |
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 | } |