xref: /linux/arch/mips/math-emu/sp_sqrt.c (revision 75bf465f0bc33e9b776a46d6a1b9b990f5fb7c37)
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 
ieee754sp_sqrt(union ieee754sp x)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 }
104