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 }