Linux Kernel  3.7.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
dfdiv.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/dfdiv.c $Revision: 1.1 $
26  *
27  * Purpose:
28  * Double Precision Floating-point Divide
29  *
30  * External Interfaces:
31  * dbl_fdiv(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 #include "float.h"
43 #include "dbl_float.h"
44 
45 /*
46  * Double Precision Floating-point Divide
47  */
48 
49 int
50 dbl_fdiv (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2,
51  dbl_floating_point * dstptr, unsigned int *status)
52 {
53  register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
54  register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
55  register int dest_exponent, count;
56  register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
57  boolean is_tiny;
58 
59  Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
60  Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
61  /*
62  * set sign bit of result
63  */
64  if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
65  Dbl_setnegativezerop1(resultp1);
66  else Dbl_setzerop1(resultp1);
67  /*
68  * check first operand for NaN's or infinity
69  */
70  if (Dbl_isinfinity_exponent(opnd1p1)) {
71  if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
72  if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
73  if (Dbl_isinfinity(opnd2p1,opnd2p2)) {
74  /*
75  * invalid since both operands
76  * are infinity
77  */
78  if (Is_invalidtrap_enabled())
79  return(INVALIDEXCEPTION);
80  Set_invalidflag();
81  Dbl_makequietnan(resultp1,resultp2);
82  Dbl_copytoptr(resultp1,resultp2,dstptr);
83  return(NOEXCEPTION);
84  }
85  /*
86  * return infinity
87  */
88  Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
89  Dbl_copytoptr(resultp1,resultp2,dstptr);
90  return(NOEXCEPTION);
91  }
92  }
93  else {
94  /*
95  * is NaN; signaling or quiet?
96  */
97  if (Dbl_isone_signaling(opnd1p1)) {
98  /* trap if INVALIDTRAP enabled */
99  if (Is_invalidtrap_enabled())
100  return(INVALIDEXCEPTION);
101  /* make NaN quiet */
102  Set_invalidflag();
103  Dbl_set_quiet(opnd1p1);
104  }
105  /*
106  * is second operand a signaling NaN?
107  */
108  else if (Dbl_is_signalingnan(opnd2p1)) {
109  /* trap if INVALIDTRAP enabled */
110  if (Is_invalidtrap_enabled())
111  return(INVALIDEXCEPTION);
112  /* make NaN quiet */
113  Set_invalidflag();
114  Dbl_set_quiet(opnd2p1);
115  Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
116  return(NOEXCEPTION);
117  }
118  /*
119  * return quiet NaN
120  */
121  Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
122  return(NOEXCEPTION);
123  }
124  }
125  /*
126  * check second operand for NaN's or infinity
127  */
128  if (Dbl_isinfinity_exponent(opnd2p1)) {
129  if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
130  /*
131  * return zero
132  */
133  Dbl_setzero_exponentmantissa(resultp1,resultp2);
134  Dbl_copytoptr(resultp1,resultp2,dstptr);
135  return(NOEXCEPTION);
136  }
137  /*
138  * is NaN; signaling or quiet?
139  */
140  if (Dbl_isone_signaling(opnd2p1)) {
141  /* trap if INVALIDTRAP enabled */
142  if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
143  /* make NaN quiet */
144  Set_invalidflag();
145  Dbl_set_quiet(opnd2p1);
146  }
147  /*
148  * return quiet NaN
149  */
150  Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
151  return(NOEXCEPTION);
152  }
153  /*
154  * check for division by zero
155  */
156  if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
157  if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
158  /* invalid since both operands are zero */
159  if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
160  Set_invalidflag();
161  Dbl_makequietnan(resultp1,resultp2);
162  Dbl_copytoptr(resultp1,resultp2,dstptr);
163  return(NOEXCEPTION);
164  }
165  if (Is_divisionbyzerotrap_enabled())
166  return(DIVISIONBYZEROEXCEPTION);
167  Set_divisionbyzeroflag();
168  Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
169  Dbl_copytoptr(resultp1,resultp2,dstptr);
170  return(NOEXCEPTION);
171  }
172  /*
173  * Generate exponent
174  */
175  dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS;
176 
177  /*
178  * Generate mantissa
179  */
180  if (Dbl_isnotzero_exponent(opnd1p1)) {
181  /* set hidden bit */
182  Dbl_clear_signexponent_set_hidden(opnd1p1);
183  }
184  else {
185  /* check for zero */
186  if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
187  Dbl_setzero_exponentmantissa(resultp1,resultp2);
188  Dbl_copytoptr(resultp1,resultp2,dstptr);
189  return(NOEXCEPTION);
190  }
191  /* is denormalized, want to normalize */
192  Dbl_clear_signexponent(opnd1p1);
193  Dbl_leftshiftby1(opnd1p1,opnd1p2);
194  Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
195  }
196  /* opnd2 needs to have hidden bit set with msb in hidden bit */
197  if (Dbl_isnotzero_exponent(opnd2p1)) {
198  Dbl_clear_signexponent_set_hidden(opnd2p1);
199  }
200  else {
201  /* is denormalized; want to normalize */
202  Dbl_clear_signexponent(opnd2p1);
203  Dbl_leftshiftby1(opnd2p1,opnd2p2);
204  while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) {
205  dest_exponent+=8;
206  Dbl_leftshiftby8(opnd2p1,opnd2p2);
207  }
208  if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) {
209  dest_exponent+=4;
210  Dbl_leftshiftby4(opnd2p1,opnd2p2);
211  }
212  while (Dbl_iszero_hidden(opnd2p1)) {
213  dest_exponent++;
214  Dbl_leftshiftby1(opnd2p1,opnd2p2);
215  }
216  }
217 
218  /* Divide the source mantissas */
219 
220  /*
221  * A non-restoring divide algorithm is used.
222  */
223  Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
224  Dbl_setzero(opnd3p1,opnd3p2);
225  for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) {
226  Dbl_leftshiftby1(opnd1p1,opnd1p2);
227  Dbl_leftshiftby1(opnd3p1,opnd3p2);
228  if (Dbl_iszero_sign(opnd1p1)) {
229  Dbl_setone_lowmantissap2(opnd3p2);
230  Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
231  }
232  else {
233  Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2);
234  }
235  }
236  if (count <= DBL_P) {
237  Dbl_leftshiftby1(opnd3p1,opnd3p2);
238  Dbl_setone_lowmantissap2(opnd3p2);
239  Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count));
240  if (Dbl_iszero_hidden(opnd3p1)) {
241  Dbl_leftshiftby1(opnd3p1,opnd3p2);
242  dest_exponent--;
243  }
244  }
245  else {
246  if (Dbl_iszero_hidden(opnd3p1)) {
247  /* need to get one more bit of result */
248  Dbl_leftshiftby1(opnd1p1,opnd1p2);
249  Dbl_leftshiftby1(opnd3p1,opnd3p2);
250  if (Dbl_iszero_sign(opnd1p1)) {
251  Dbl_setone_lowmantissap2(opnd3p2);
252  Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
253  }
254  else {
255  Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
256  }
257  dest_exponent--;
258  }
259  if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE;
260  stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2);
261  }
262  inexact = guardbit | stickybit;
263 
264  /*
265  * round result
266  */
267  if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
268  Dbl_clear_signexponent(opnd3p1);
269  switch (Rounding_mode()) {
270  case ROUNDPLUS:
271  if (Dbl_iszero_sign(resultp1))
272  Dbl_increment(opnd3p1,opnd3p2);
273  break;
274  case ROUNDMINUS:
275  if (Dbl_isone_sign(resultp1))
276  Dbl_increment(opnd3p1,opnd3p2);
277  break;
278  case ROUNDNEAREST:
279  if (guardbit && (stickybit ||
280  Dbl_isone_lowmantissap2(opnd3p2))) {
281  Dbl_increment(opnd3p1,opnd3p2);
282  }
283  }
284  if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
285  }
286  Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);
287 
288  /*
289  * Test for overflow
290  */
291  if (dest_exponent >= DBL_INFINITY_EXPONENT) {
292  /* trap if OVERFLOWTRAP enabled */
293  if (Is_overflowtrap_enabled()) {
294  /*
295  * Adjust bias of result
296  */
297  Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
298  Dbl_copytoptr(resultp1,resultp2,dstptr);
299  if (inexact)
300  if (Is_inexacttrap_enabled())
301  return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
302  else Set_inexactflag();
303  return(OVERFLOWEXCEPTION);
304  }
305  Set_overflowflag();
306  /* set result to infinity or largest number */
307  Dbl_setoverflow(resultp1,resultp2);
308  inexact = TRUE;
309  }
310  /*
311  * Test for underflow
312  */
313  else if (dest_exponent <= 0) {
314  /* trap if UNDERFLOWTRAP enabled */
315  if (Is_underflowtrap_enabled()) {
316  /*
317  * Adjust bias of result
318  */
319  Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
320  Dbl_copytoptr(resultp1,resultp2,dstptr);
321  if (inexact)
322  if (Is_inexacttrap_enabled())
323  return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
324  else Set_inexactflag();
325  return(UNDERFLOWEXCEPTION);
326  }
327 
328  /* Determine if should set underflow flag */
329  is_tiny = TRUE;
330  if (dest_exponent == 0 && inexact) {
331  switch (Rounding_mode()) {
332  case ROUNDPLUS:
333  if (Dbl_iszero_sign(resultp1)) {
334  Dbl_increment(opnd3p1,opnd3p2);
335  if (Dbl_isone_hiddenoverflow(opnd3p1))
336  is_tiny = FALSE;
337  Dbl_decrement(opnd3p1,opnd3p2);
338  }
339  break;
340  case ROUNDMINUS:
341  if (Dbl_isone_sign(resultp1)) {
342  Dbl_increment(opnd3p1,opnd3p2);
343  if (Dbl_isone_hiddenoverflow(opnd3p1))
344  is_tiny = FALSE;
345  Dbl_decrement(opnd3p1,opnd3p2);
346  }
347  break;
348  case ROUNDNEAREST:
349  if (guardbit && (stickybit ||
350  Dbl_isone_lowmantissap2(opnd3p2))) {
351  Dbl_increment(opnd3p1,opnd3p2);
352  if (Dbl_isone_hiddenoverflow(opnd3p1))
353  is_tiny = FALSE;
354  Dbl_decrement(opnd3p1,opnd3p2);
355  }
356  break;
357  }
358  }
359 
360  /*
361  * denormalize result or set to signed zero
362  */
363  stickybit = inexact;
364  Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
365  stickybit,inexact);
366 
367  /* return rounded number */
368  if (inexact) {
369  switch (Rounding_mode()) {
370  case ROUNDPLUS:
371  if (Dbl_iszero_sign(resultp1)) {
372  Dbl_increment(opnd3p1,opnd3p2);
373  }
374  break;
375  case ROUNDMINUS:
376  if (Dbl_isone_sign(resultp1)) {
377  Dbl_increment(opnd3p1,opnd3p2);
378  }
379  break;
380  case ROUNDNEAREST:
381  if (guardbit && (stickybit ||
382  Dbl_isone_lowmantissap2(opnd3p2))) {
383  Dbl_increment(opnd3p1,opnd3p2);
384  }
385  break;
386  }
387  if (is_tiny) Set_underflowflag();
388  }
389  Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
390  }
391  else Dbl_set_exponent(resultp1,dest_exponent);
392  Dbl_copytoptr(resultp1,resultp2,dstptr);
393 
394  /* check for inexact */
395  if (inexact) {
396  if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
397  else Set_inexactflag();
398  }
399  return(NOEXCEPTION);
400 }
Generated on Thu Jan 10 2013 13:12:59 for Linux Kernel by   doxygen 1.8.2