/arch/parisc/math-emu/dfrem.c
C | 297 lines | 149 code | 14 blank | 134 comment | 47 complexity | 3f748efc3c2abe1e256169e91e6a0d82 MD5 | raw file
Possible License(s): LGPL-2.0, AGPL-1.0, GPL-2.0
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
50int
51dbl_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}