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
12 union 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 }