Linux Kernel  3.7.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
dfrem.c
Go to the documentation of this file.
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) <[email protected]>
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 
50 int
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 */
69  return(INVALIDEXCEPTION);
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 */
83  return(INVALIDEXCEPTION);
84  /* make NaN quiet */
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 */
94  return(INVALIDEXCEPTION);
95  /* make NaN quiet */
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 */
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 */
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 {
167  }
168  if (opnd2_exponent == 0) {
169  /* normalize, then continue */
170  opnd2_exponent = 1;
171  Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
172  }
173  else {
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 }