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