File: /usr/src/linux/arch/mips/math-emu/sp_sqrt.c

1     /* IEEE754 floating point arithmetic
2      * single precision square root
3      */
4     /*
5      * MIPS floating point support
6      * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
7      * http://www.algor.co.uk
8      *
9      * ########################################################################
10      *
11      *  This program is free software; you can distribute it and/or modify it
12      *  under the terms of the GNU General Public License (Version 2) as
13      *  published by the Free Software Foundation.
14      *
15      *  This program is distributed in the hope it will be useful, but WITHOUT
16      *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17      *  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
18      *  for more details.
19      *
20      *  You should have received a copy of the GNU General Public License along
21      *  with this program; if not, write to the Free Software Foundation, Inc.,
22      *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23      *
24      * ########################################################################
25      */
26     
27     
28     #include "ieee754sp.h"
29     
30     static const struct ieee754sp_konst knan = {
31     	0, SP_EBIAS + SP_EMAX + 1, 0
32     };
33     
34     #define nan	((ieee754sp)knan)
35     
36     ieee754sp ieee754sp_sqrt(ieee754sp x)
37     {
38     	int sign = (int) 0x80000000;
39     	int ix, s, q, m, t, i;
40     	unsigned int r;
41     	COMPXDP;
42     
43     	/* take care of Inf and NaN */
44     
45     	EXPLODEXDP;
46     
47     	/* x == INF or NAN? */
48     	switch (xc) {
49     	case IEEE754_CLASS_QNAN:
50     	case IEEE754_CLASS_SNAN:
51     		/* sqrt(Nan) = Nan */
52     		return ieee754sp_nanxcpt(x, "sqrt");
53     	case IEEE754_CLASS_ZERO:
54     		/* sqrt(0) = 0 */
55     		return x;
56     	case IEEE754_CLASS_INF:
57     		if (xs)
58     			/* sqrt(-Inf) = Nan */
59     			return ieee754sp_nanxcpt(nan, "sqrt");
60     		/* sqrt(+Inf) = Inf */
61     		return x;
62     	case IEEE754_CLASS_DNORM:
63     	case IEEE754_CLASS_NORM:
64     		if (xs)
65     			/* sqrt(-x) = Nan */
66     			return ieee754sp_nanxcpt(nan, "sqrt");
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     		switch (ieee754_csr.rm) {
103     		case IEEE754_RP:
104     			q += 2;
105     			break;
106     		case IEEE754_RN:
107     			q += (q & 1);
108     			break;
109     		}
110     	}
111     	ix = (q >> 1) + 0x3f000000;
112     	ix += (m << 23);
113     	x.bits = ix;
114     	return x;
115     }
116