Linux Kernel  3.7.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
sp_sqrt.c
Go to the documentation of this file.
1 /* IEEE754 floating point arithmetic
2  * single precision square root
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 "ieee754sp.h"
28 
29 ieee754sp ieee754sp_sqrt(ieee754sp x)
30 {
31  int ix, s, q, m, t, i;
32  unsigned int r;
33  COMPXSP;
34 
35  /* take care of Inf and NaN */
36 
37  EXPLODEXSP;
38  CLEARCX;
39  FLUSHXSP;
40 
41  /* x == INF or NAN? */
42  switch (xc) {
43  case IEEE754_CLASS_QNAN:
44  /* sqrt(Nan) = Nan */
45  return ieee754sp_nanxcpt(x, "sqrt");
46  case IEEE754_CLASS_SNAN:
47  SETCX(IEEE754_INVALID_OPERATION);
48  return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
49  case IEEE754_CLASS_ZERO:
50  /* sqrt(0) = 0 */
51  return x;
52  case IEEE754_CLASS_INF:
53  if (xs) {
54  /* sqrt(-Inf) = Nan */
55  SETCX(IEEE754_INVALID_OPERATION);
56  return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
57  }
58  /* sqrt(+Inf) = Inf */
59  return x;
60  case IEEE754_CLASS_DNORM:
61  case IEEE754_CLASS_NORM:
62  if (xs) {
63  /* sqrt(-x) = Nan */
64  SETCX(IEEE754_INVALID_OPERATION);
65  return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
66  }
67  break;
68  }
69 
70  ix = x.bits;
71 
72  /* normalize x */
73  m = (ix >> 23);
74  if (m == 0) { /* subnormal x */
75  for (i = 0; (ix & 0x00800000) == 0; i++)
76  ix <<= 1;
77  m -= i - 1;
78  }
79  m -= 127; /* unbias exponent */
80  ix = (ix & 0x007fffff) | 0x00800000;
81  if (m & 1) /* odd m, double x to make it even */
82  ix += ix;
83  m >>= 1; /* m = [m/2] */
84 
85  /* generate sqrt(x) bit by bit */
86  ix += ix;
87  q = s = 0; /* q = sqrt(x) */
88  r = 0x01000000; /* r = moving bit from right to left */
89 
90  while (r != 0) {
91  t = s + r;
92  if (t <= ix) {
93  s = t + r;
94  ix -= t;
95  q += r;
96  }
97  ix += ix;
98  r >>= 1;
99  }
100 
101  if (ix != 0) {
102  SETCX(IEEE754_INEXACT);
103  switch (ieee754_csr.rm) {
104  case IEEE754_RP:
105  q += 2;
106  break;
107  case IEEE754_RN:
108  q += (q & 1);
109  break;
110  }
111  }
112  ix = (q >> 1) + 0x3f000000;
113  ix += (m << 23);
114  x.bits = ix;
115  return x;
116 }
Generated on Thu Jan 10 2013 13:11:40 for Linux Kernel by   doxygen 1.8.2