sp_sqrt.c

  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}