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}
104