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// SPDX-License-Identifier: GPL-2.0-only
2/* IEEE754 floating point arithmetic
3 * single precision square root
4 */
5/*
6 * MIPS floating point support
7 * Copyright (C) 1994-2000 Algorithmics Ltd.
8 */
9
10#include "ieee754sp.h"
11
12union ieee754sp ieee754sp_sqrt(union ieee754sp x)
13{
14 int ix, s, q, m, t, i;
15 unsigned int r;
16 COMPXSP;
17
18 /* take care of Inf and NaN */
19
20 EXPLODEXSP;
21 ieee754_clearcx();
22 FLUSHXSP;
23
24 /* x == INF or NAN? */
25 switch (xc) {
26 case IEEE754_CLASS_SNAN:
27 return ieee754sp_nanxcpt(x);
28
29 case IEEE754_CLASS_QNAN:
30 /* sqrt(Nan) = Nan */
31 return x;
32
33 case IEEE754_CLASS_ZERO:
34 /* sqrt(0) = 0 */
35 return x;
36
37 case IEEE754_CLASS_INF:
38 if (xs) {
39 /* sqrt(-Inf) = Nan */
40 ieee754_setcx(IEEE754_INVALID_OPERATION);
41 return ieee754sp_indef();
42 }
43 /* sqrt(+Inf) = Inf */
44 return x;
45
46 case IEEE754_CLASS_DNORM:
47 case IEEE754_CLASS_NORM:
48 if (xs) {
49 /* sqrt(-x) = Nan */
50 ieee754_setcx(IEEE754_INVALID_OPERATION);
51 return ieee754sp_indef();
52 }
53 break;
54 }
55
56 ix = x.bits;
57
58 /* normalize x */
59 m = (ix >> 23);
60 if (m == 0) { /* subnormal x */
61 for (i = 0; (ix & 0x00800000) == 0; i++)
62 ix <<= 1;
63 m -= i - 1;
64 }
65 m -= 127; /* unbias exponent */
66 ix = (ix & 0x007fffff) | 0x00800000;
67 if (m & 1) /* odd m, double x to make it even */
68 ix += ix;
69 m >>= 1; /* m = [m/2] */
70
71 /* generate sqrt(x) bit by bit */
72 ix += ix;
73 s = 0;
74 q = 0; /* q = sqrt(x) */
75 r = 0x01000000; /* r = moving bit from right to left */
76
77 while (r != 0) {
78 t = s + r;
79 if (t <= ix) {
80 s = t + r;
81 ix -= t;
82 q += r;
83 }
84 ix += ix;
85 r >>= 1;
86 }
87
88 if (ix != 0) {
89 ieee754_setcx(IEEE754_INEXACT);
90 switch (ieee754_csr.rm) {
91 case FPU_CSR_RU:
92 q += 2;
93 break;
94 case FPU_CSR_RN:
95 q += (q & 1);
96 break;
97 }
98 }
99 ix = (q >> 1) + 0x3f000000;
100 ix += (m << 23);
101 x.bits = ix;
102 return x;
103}