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., 675 Mass Ave, Cambridge, MA 02139, USA.
20 */
21 /*
22 * BEGIN_DESC
23 *
24 * File:
25 * @(#) pa/spmath/sfrem.c $Revision: $
26 *
27 * Purpose:
28 * Single Precision Floating-point Remainder
29 *
30 * External Interfaces:
31 * sgl_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 "sgl_float.h"
45
46 /*
47 * Single Precision Floating-point Remainder
48 */
49
50 int
sgl_frem(sgl_floating_point * srcptr1,sgl_floating_point * srcptr2,sgl_floating_point * dstptr,unsigned int * status)51 sgl_frem (sgl_floating_point * srcptr1, sgl_floating_point * srcptr2,
52 sgl_floating_point * dstptr, unsigned int *status)
53 {
54 register unsigned int opnd1, opnd2, result;
55 register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
56 register boolean roundup = FALSE;
57
58 opnd1 = *srcptr1;
59 opnd2 = *srcptr2;
60 /*
61 * check first operand for NaN's or infinity
62 */
63 if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) {
64 if (Sgl_iszero_mantissa(opnd1)) {
65 if (Sgl_isnotnan(opnd2)) {
66 /* invalid since first operand is infinity */
67 if (Is_invalidtrap_enabled())
68 return(INVALIDEXCEPTION);
69 Set_invalidflag();
70 Sgl_makequietnan(result);
71 *dstptr = result;
72 return(NOEXCEPTION);
73 }
74 }
75 else {
76 /*
77 * is NaN; signaling or quiet?
78 */
79 if (Sgl_isone_signaling(opnd1)) {
80 /* trap if INVALIDTRAP enabled */
81 if (Is_invalidtrap_enabled())
82 return(INVALIDEXCEPTION);
83 /* make NaN quiet */
84 Set_invalidflag();
85 Sgl_set_quiet(opnd1);
86 }
87 /*
88 * is second operand a signaling NaN?
89 */
90 else if (Sgl_is_signalingnan(opnd2)) {
91 /* trap if INVALIDTRAP enabled */
92 if (Is_invalidtrap_enabled())
93 return(INVALIDEXCEPTION);
94 /* make NaN quiet */
95 Set_invalidflag();
96 Sgl_set_quiet(opnd2);
97 *dstptr = opnd2;
98 return(NOEXCEPTION);
99 }
100 /*
101 * return quiet NaN
102 */
103 *dstptr = opnd1;
104 return(NOEXCEPTION);
105 }
106 }
107 /*
108 * check second operand for NaN's or infinity
109 */
110 if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) {
111 if (Sgl_iszero_mantissa(opnd2)) {
112 /*
113 * return first operand
114 */
115 *dstptr = opnd1;
116 return(NOEXCEPTION);
117 }
118 /*
119 * is NaN; signaling or quiet?
120 */
121 if (Sgl_isone_signaling(opnd2)) {
122 /* trap if INVALIDTRAP enabled */
123 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
124 /* make NaN quiet */
125 Set_invalidflag();
126 Sgl_set_quiet(opnd2);
127 }
128 /*
129 * return quiet NaN
130 */
131 *dstptr = opnd2;
132 return(NOEXCEPTION);
133 }
134 /*
135 * check second operand for zero
136 */
137 if (Sgl_iszero_exponentmantissa(opnd2)) {
138 /* invalid since second operand is zero */
139 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
140 Set_invalidflag();
141 Sgl_makequietnan(result);
142 *dstptr = result;
143 return(NOEXCEPTION);
144 }
145
146 /*
147 * get sign of result
148 */
149 result = opnd1;
150
151 /*
152 * check for denormalized operands
153 */
154 if (opnd1_exponent == 0) {
155 /* check for zero */
156 if (Sgl_iszero_mantissa(opnd1)) {
157 *dstptr = opnd1;
158 return(NOEXCEPTION);
159 }
160 /* normalize, then continue */
161 opnd1_exponent = 1;
162 Sgl_normalize(opnd1,opnd1_exponent);
163 }
164 else {
165 Sgl_clear_signexponent_set_hidden(opnd1);
166 }
167 if (opnd2_exponent == 0) {
168 /* normalize, then continue */
169 opnd2_exponent = 1;
170 Sgl_normalize(opnd2,opnd2_exponent);
171 }
172 else {
173 Sgl_clear_signexponent_set_hidden(opnd2);
174 }
175
176 /* find result exponent and divide step loop count */
177 dest_exponent = opnd2_exponent - 1;
178 stepcount = opnd1_exponent - opnd2_exponent;
179
180 /*
181 * check for opnd1/opnd2 < 1
182 */
183 if (stepcount < 0) {
184 /*
185 * check for opnd1/opnd2 > 1/2
186 *
187 * In this case n will round to 1, so
188 * r = opnd1 - opnd2
189 */
190 if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) {
191 Sgl_all(result) = ~Sgl_all(result); /* set sign */
192 /* align opnd2 with opnd1 */
193 Sgl_leftshiftby1(opnd2);
194 Sgl_subtract(opnd2,opnd1,opnd2);
195 /* now normalize */
196 while (Sgl_iszero_hidden(opnd2)) {
197 Sgl_leftshiftby1(opnd2);
198 dest_exponent--;
199 }
200 Sgl_set_exponentmantissa(result,opnd2);
201 goto testforunderflow;
202 }
203 /*
204 * opnd1/opnd2 <= 1/2
205 *
206 * In this case n will round to zero, so
207 * r = opnd1
208 */
209 Sgl_set_exponentmantissa(result,opnd1);
210 dest_exponent = opnd1_exponent;
211 goto testforunderflow;
212 }
213
214 /*
215 * Generate result
216 *
217 * Do iterative subtract until remainder is less than operand 2.
218 */
219 while (stepcount-- > 0 && Sgl_all(opnd1)) {
220 if (Sgl_isnotlessthan(opnd1,opnd2))
221 Sgl_subtract(opnd1,opnd2,opnd1);
222 Sgl_leftshiftby1(opnd1);
223 }
224 /*
225 * Do last subtract, then determine which way to round if remainder
226 * is exactly 1/2 of opnd2
227 */
228 if (Sgl_isnotlessthan(opnd1,opnd2)) {
229 Sgl_subtract(opnd1,opnd2,opnd1);
230 roundup = TRUE;
231 }
232 if (stepcount > 0 || Sgl_iszero(opnd1)) {
233 /* division is exact, remainder is zero */
234 Sgl_setzero_exponentmantissa(result);
235 *dstptr = result;
236 return(NOEXCEPTION);
237 }
238
239 /*
240 * Check for cases where opnd1/opnd2 < n
241 *
242 * In this case the result's sign will be opposite that of
243 * opnd1. The mantissa also needs some correction.
244 */
245 Sgl_leftshiftby1(opnd1);
246 if (Sgl_isgreaterthan(opnd1,opnd2)) {
247 Sgl_invert_sign(result);
248 Sgl_subtract((opnd2<<1),opnd1,opnd1);
249 }
250 /* check for remainder being exactly 1/2 of opnd2 */
251 else if (Sgl_isequal(opnd1,opnd2) && roundup) {
252 Sgl_invert_sign(result);
253 }
254
255 /* normalize result's mantissa */
256 while (Sgl_iszero_hidden(opnd1)) {
257 dest_exponent--;
258 Sgl_leftshiftby1(opnd1);
259 }
260 Sgl_set_exponentmantissa(result,opnd1);
261
262 /*
263 * Test for underflow
264 */
265 testforunderflow:
266 if (dest_exponent <= 0) {
267 /* trap if UNDERFLOWTRAP enabled */
268 if (Is_underflowtrap_enabled()) {
269 /*
270 * Adjust bias of result
271 */
272 Sgl_setwrapped_exponent(result,dest_exponent,unfl);
273 *dstptr = result;
274 /* frem is always exact */
275 return(UNDERFLOWEXCEPTION);
276 }
277 /*
278 * denormalize result or set to signed zero
279 */
280 if (dest_exponent >= (1 - SGL_P)) {
281 Sgl_rightshift_exponentmantissa(result,1-dest_exponent);
282 }
283 else {
284 Sgl_setzero_exponentmantissa(result);
285 }
286 }
287 else Sgl_set_exponent(result,dest_exponent);
288 *dstptr = result;
289 return(NOEXCEPTION);
290 }
291