Loading...
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
29ieee754sp 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}
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
29ieee754sp 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}