Linux Kernel  3.7.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
dp_sub.c
Go to the documentation of this file.
1 /* IEEE754 floating point arithmetic
2  * double precision: common utilities
3  */
4 /*
5  * MIPS floating point support
6  * Copyright (C) 1994-2000 Algorithmics Ltd.
7  *
8  * ########################################################################
9  *
10  * This program is free software; you can distribute it and/or modify it
11  * under the terms of the GNU General Public License (Version 2) as
12  * published by the Free Software Foundation.
13  *
14  * This program is distributed in the hope it will be useful, but WITHOUT
15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17  * for more details.
18  *
19  * You should have received a copy of the GNU General Public License along
20  * with this program; if not, write to the Free Software Foundation, Inc.,
21  * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
22  *
23  * ########################################################################
24  */
25 
26 
27 #include "ieee754dp.h"
28 
29 ieee754dp ieee754dp_sub(ieee754dp x, ieee754dp y)
30 {
31  COMPXDP;
32  COMPYDP;
33 
34  EXPLODEXDP;
35  EXPLODEYDP;
36 
37  CLEARCX;
38 
39  FLUSHXDP;
40  FLUSHYDP;
41 
42  switch (CLPAIR(xc, yc)) {
55  return ieee754dp_nanxcpt(ieee754dp_indef(), "sub", x, y);
56 
61  return y;
62 
68  return x;
69 
70 
71  /* Infinity handling
72  */
73 
75  if (xs != ys)
76  return x;
78  return ieee754dp_xcpt(ieee754dp_indef(), "sub", x, y);
79 
83  return ieee754dp_inf(ys ^ 1);
84 
88  return x;
89 
90  /* Zero handling
91  */
92 
94  if (xs != ys)
95  return x;
96  else
97  return ieee754dp_zero(ieee754_csr.rm ==
98  IEEE754_RD);
99 
102  return x;
103 
106  /* quick fix up */
107  DPSIGN(y) ^= 1;
108  return y;
109 
111  DPDNORMX;
112  /* FALL THROUGH */
113 
115  /* normalize ym,ye */
116  DPDNORMY;
117  break;
118 
120  /* normalize xm,xe */
121  DPDNORMX;
122  break;
123 
125  break;
126  }
127  /* flip sign of y and handle as add */
128  ys ^= 1;
129 
130  assert(xm & DP_HIDDEN_BIT);
131  assert(ym & DP_HIDDEN_BIT);
132 
133 
134  /* provide guard,round and stick bit dpace */
135  xm <<= 3;
136  ym <<= 3;
137 
138  if (xe > ye) {
139  /* have to shift y fraction right to align
140  */
141  int s = xe - ye;
142  ym = XDPSRS(ym, s);
143  ye += s;
144  } else if (ye > xe) {
145  /* have to shift x fraction right to align
146  */
147  int s = ye - xe;
148  xm = XDPSRS(xm, s);
149  xe += s;
150  }
151  assert(xe == ye);
152  assert(xe <= DP_EMAX);
153 
154  if (xs == ys) {
155  /* generate 28 bit result of adding two 27 bit numbers
156  */
157  xm = xm + ym;
158  xe = xe;
159  xs = xs;
160 
161  if (xm >> (DP_MBITS + 1 + 3)) { /* carry out */
162  xm = XDPSRS1(xm); /* shift preserving sticky */
163  xe++;
164  }
165  } else {
166  if (xm >= ym) {
167  xm = xm - ym;
168  xe = xe;
169  xs = xs;
170  } else {
171  xm = ym - xm;
172  xe = xe;
173  xs = ys;
174  }
175  if (xm == 0) {
176  if (ieee754_csr.rm == IEEE754_RD)
177  return ieee754dp_zero(1); /* round negative inf. => sign = -1 */
178  else
179  return ieee754dp_zero(0); /* other round modes => sign = 1 */
180  }
181 
182  /* normalize to rounding precision
183  */
184  while ((xm >> (DP_MBITS + 3)) == 0) {
185  xm <<= 1;
186  xe--;
187  }
188  }
189  DPNORMRET2(xs, xe, xm, "sub", x, y);
190 }