Linux Kernel  3.7.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
sfsqrt.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/sfsqrt.c $Revision: 1.1 $
26  *
27  * Purpose:
28  * Single Floating-point Square Root
29  *
30  * External Interfaces:
31  * sgl_fsqrt(srcptr,nullptr,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 "sgl_float.h"
44 
45 /*
46  * Single Floating-point Square Root
47  */
48 
49 /*ARGSUSED*/
50 unsigned int
52  sgl_floating_point *srcptr,
53  unsigned int *nullptr,
54  sgl_floating_point *dstptr,
55  unsigned int *status)
56 {
57  register unsigned int src, result;
58  register int src_exponent;
59  register unsigned int newbit, sum;
60  register boolean guardbit = FALSE, even_exponent;
61 
62  src = *srcptr;
63  /*
64  * check source operand for NaN or infinity
65  */
66  if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) {
67  /*
68  * is signaling NaN?
69  */
70  if (Sgl_isone_signaling(src)) {
71  /* trap if INVALIDTRAP enabled */
73  /* make NaN quiet */
75  Sgl_set_quiet(src);
76  }
77  /*
78  * Return quiet NaN or positive infinity.
79  * Fall through to negative test if negative infinity.
80  */
81  if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) {
82  *dstptr = src;
83  return(NOEXCEPTION);
84  }
85  }
86 
87  /*
88  * check for zero source operand
89  */
90  if (Sgl_iszero_exponentmantissa(src)) {
91  *dstptr = src;
92  return(NOEXCEPTION);
93  }
94 
95  /*
96  * check for negative source operand
97  */
98  if (Sgl_isone_sign(src)) {
99  /* trap if INVALIDTRAP enabled */
101  /* make NaN quiet */
102  Set_invalidflag();
103  Sgl_makequietnan(src);
104  *dstptr = src;
105  return(NOEXCEPTION);
106  }
107 
108  /*
109  * Generate result
110  */
111  if (src_exponent > 0) {
112  even_exponent = Sgl_hidden(src);
114  }
115  else {
116  /* normalize operand */
118  src_exponent++;
119  Sgl_normalize(src,src_exponent);
120  even_exponent = src_exponent & 1;
121  }
122  if (even_exponent) {
123  /* exponent is even */
124  /* Add comment here. Explain why odd exponent needs correction */
125  Sgl_leftshiftby1(src);
126  }
127  /*
128  * Add comment here. Explain following algorithm.
129  *
130  * Trust me, it works.
131  *
132  */
133  Sgl_setzero(result);
134  newbit = 1 << SGL_P;
135  while (newbit && Sgl_isnotzero(src)) {
136  Sgl_addition(result,newbit,sum);
137  if(sum <= Sgl_all(src)) {
138  /* update result */
139  Sgl_addition(result,(newbit<<1),result);
140  Sgl_subtract(src,sum,src);
141  }
142  Sgl_rightshiftby1(newbit);
143  Sgl_leftshiftby1(src);
144  }
145  /* correct exponent for pre-shift */
146  if (even_exponent) {
147  Sgl_rightshiftby1(result);
148  }
149 
150  /* check for inexact */
151  if (Sgl_isnotzero(src)) {
152  if (!even_exponent && Sgl_islessthan(result,src))
153  Sgl_increment(result);
154  guardbit = Sgl_lowmantissa(result);
155  Sgl_rightshiftby1(result);
156 
157  /* now round result */
158  switch (Rounding_mode()) {
159  case ROUNDPLUS:
160  Sgl_increment(result);
161  break;
162  case ROUNDNEAREST:
163  /* stickybit is always true, so guardbit
164  * is enough to determine rounding */
165  if (guardbit) {
166  Sgl_increment(result);
167  }
168  break;
169  }
170  /* increment result exponent by 1 if mantissa overflowed */
171  if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2;
172 
173  if (Is_inexacttrap_enabled()) {
174  Sgl_set_exponent(result,
175  ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
176  *dstptr = result;
177  return(INEXACTEXCEPTION);
178  }
179  else Set_inexactflag();
180  }
181  else {
182  Sgl_rightshiftby1(result);
183  }
184  Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
185  *dstptr = result;
186  return(NOEXCEPTION);
187 }