Linux 2.6.25-rc4
[linux-2.6/next.git] / arch / parisc / math-emu / dfrem.c
blobb983785346778bda4999d67f7fe888f0cab19a76
1 /*
2 * Linux/PA-RISC Project (http://www.parisc-linux.org/)
4 * Floating-point emulation code
5 * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
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.
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.
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
22 * BEGIN_DESC
24 * File:
25 * @(#) pa/spmath/dfrem.c $Revision: 1.1 $
27 * Purpose:
28 * Double Precision Floating-point Remainder
30 * External Interfaces:
31 * dbl_frem(srcptr1,srcptr2,dstptr,status)
33 * Internal Interfaces:
35 * Theory:
36 * <<please update with a overview of the operation of this file>>
38 * END_DESC
43 #include "float.h"
44 #include "dbl_float.h"
47 * Double Precision Floating-point Remainder
50 int
51 dbl_frem (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2,
52 dbl_floating_point * dstptr, unsigned int *status)
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;
59 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
60 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
62 * check first operand for NaN's or infinity
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);
76 else {
78 * is NaN; signaling or quiet?
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);
88 /*
89 * is second operand a signaling NaN?
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);
102 * return quiet NaN
104 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
105 return(NOEXCEPTION);
109 * check second operand for NaN's or infinity
111 if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
112 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
114 * return first operand
116 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
117 return(NOEXCEPTION);
120 * is NaN; signaling or quiet?
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);
130 * return quiet NaN
132 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
133 return(NOEXCEPTION);
136 * check second operand for zero
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);
148 * get sign of result
150 resultp1 = opnd1p1;
153 * check for denormalized operands
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);
161 /* normalize, then continue */
162 opnd1_exponent = 1;
163 Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
165 else {
166 Dbl_clear_signexponent_set_hidden(opnd1p1);
168 if (opnd2_exponent == 0) {
169 /* normalize, then continue */
170 opnd2_exponent = 1;
171 Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
173 else {
174 Dbl_clear_signexponent_set_hidden(opnd2p1);
177 /* find result exponent and divide step loop count */
178 dest_exponent = opnd2_exponent - 1;
179 stepcount = opnd1_exponent - opnd2_exponent;
182 * check for opnd1/opnd2 < 1
184 if (stepcount < 0) {
186 * check for opnd1/opnd2 > 1/2
188 * In this case n will round to 1, so
189 * r = opnd1 - opnd2
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--;
204 Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
205 goto testforunderflow;
208 * opnd1/opnd2 <= 1/2
210 * In this case n will round to zero, so
211 * r = opnd1
213 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
214 dest_exponent = opnd1_exponent;
215 goto testforunderflow;
219 * Generate result
221 * Do iterative subtract until remainder is less than operand 2.
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);
227 Dbl_leftshiftby1(opnd1p1,opnd1p2);
230 * Do last subtract, then determine which way to round if remainder
231 * is exactly 1/2 of opnd2
233 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
234 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
235 roundup = TRUE;
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);
245 * Check for cases where opnd1/opnd2 < n
247 * In this case the result's sign will be opposite that of
248 * opnd1. The mantissa also needs some correction.
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);
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);
261 /* normalize result's mantissa */
262 while (Dbl_iszero_hidden(opnd1p1)) {
263 dest_exponent--;
264 Dbl_leftshiftby1(opnd1p1,opnd1p2);
266 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
269 * Test for underflow
271 testforunderflow:
272 if (dest_exponent <= 0) {
273 /* trap if UNDERFLOWTRAP enabled */
274 if (Is_underflowtrap_enabled()) {
276 * Adjust bias of result
278 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
279 /* frem is always exact */
280 Dbl_copytoptr(resultp1,resultp2,dstptr);
281 return(UNDERFLOWEXCEPTION);
284 * denormalize result or set to signed zero
286 if (dest_exponent >= (1 - DBL_P)) {
287 Dbl_rightshift_exponentmantissa(resultp1,resultp2,
288 1-dest_exponent);
290 else {
291 Dbl_setzero_exponentmantissa(resultp1,resultp2);
294 else Dbl_set_exponent(resultp1,dest_exponent);
295 Dbl_copytoptr(resultp1,resultp2,dstptr);
296 return(NOEXCEPTION);