xref: /openbmc/linux/arch/mips/math-emu/sp_sqrt.c (revision 90efba36)
11da177e4SLinus Torvalds /* IEEE754 floating point arithmetic
21da177e4SLinus Torvalds  * single precision square root
31da177e4SLinus Torvalds  */
41da177e4SLinus Torvalds /*
51da177e4SLinus Torvalds  * MIPS floating point support
61da177e4SLinus Torvalds  * Copyright (C) 1994-2000 Algorithmics Ltd.
71da177e4SLinus Torvalds  *
81da177e4SLinus Torvalds  * ########################################################################
91da177e4SLinus Torvalds  *
101da177e4SLinus Torvalds  *  This program is free software; you can distribute it and/or modify it
111da177e4SLinus Torvalds  *  under the terms of the GNU General Public License (Version 2) as
121da177e4SLinus Torvalds  *  published by the Free Software Foundation.
131da177e4SLinus Torvalds  *
141da177e4SLinus Torvalds  *  This program is distributed in the hope it will be useful, but WITHOUT
151da177e4SLinus Torvalds  *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
161da177e4SLinus Torvalds  *  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
171da177e4SLinus Torvalds  *  for more details.
181da177e4SLinus Torvalds  *
191da177e4SLinus Torvalds  *  You should have received a copy of the GNU General Public License along
201da177e4SLinus Torvalds  *  with this program; if not, write to the Free Software Foundation, Inc.,
211da177e4SLinus Torvalds  *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
221da177e4SLinus Torvalds  *
231da177e4SLinus Torvalds  * ########################################################################
241da177e4SLinus Torvalds  */
251da177e4SLinus Torvalds 
261da177e4SLinus Torvalds 
271da177e4SLinus Torvalds #include "ieee754sp.h"
281da177e4SLinus Torvalds 
292209bcb1SRalf Baechle union ieee754sp ieee754sp_sqrt(union ieee754sp x)
301da177e4SLinus Torvalds {
311da177e4SLinus Torvalds 	int ix, s, q, m, t, i;
321da177e4SLinus Torvalds 	unsigned int r;
331da177e4SLinus Torvalds 	COMPXSP;
341da177e4SLinus Torvalds 
351da177e4SLinus Torvalds 	/* take care of Inf and NaN */
361da177e4SLinus Torvalds 
371da177e4SLinus Torvalds 	EXPLODEXSP;
389e8bad1fSRalf Baechle 	ieee754_clearcx();
391da177e4SLinus Torvalds 	FLUSHXSP;
401da177e4SLinus Torvalds 
411da177e4SLinus Torvalds 	/* x == INF or NAN? */
421da177e4SLinus Torvalds 	switch (xc) {
431da177e4SLinus Torvalds 	case IEEE754_CLASS_QNAN:
441da177e4SLinus Torvalds 		/* sqrt(Nan) = Nan */
4590efba36SRalf Baechle 		return ieee754sp_nanxcpt(x);
461da177e4SLinus Torvalds 	case IEEE754_CLASS_SNAN:
479e8bad1fSRalf Baechle 		ieee754_setcx(IEEE754_INVALID_OPERATION);
4890efba36SRalf Baechle 		return ieee754sp_nanxcpt(ieee754sp_indef());
491da177e4SLinus Torvalds 	case IEEE754_CLASS_ZERO:
501da177e4SLinus Torvalds 		/* sqrt(0) = 0 */
511da177e4SLinus Torvalds 		return x;
521da177e4SLinus Torvalds 	case IEEE754_CLASS_INF:
531da177e4SLinus Torvalds 		if (xs) {
541da177e4SLinus Torvalds 			/* sqrt(-Inf) = Nan */
559e8bad1fSRalf Baechle 			ieee754_setcx(IEEE754_INVALID_OPERATION);
5690efba36SRalf Baechle 			return ieee754sp_nanxcpt(ieee754sp_indef());
571da177e4SLinus Torvalds 		}
581da177e4SLinus Torvalds 		/* sqrt(+Inf) = Inf */
591da177e4SLinus Torvalds 		return x;
601da177e4SLinus Torvalds 	case IEEE754_CLASS_DNORM:
611da177e4SLinus Torvalds 	case IEEE754_CLASS_NORM:
621da177e4SLinus Torvalds 		if (xs) {
631da177e4SLinus Torvalds 			/* sqrt(-x) = Nan */
649e8bad1fSRalf Baechle 			ieee754_setcx(IEEE754_INVALID_OPERATION);
6590efba36SRalf Baechle 			return ieee754sp_nanxcpt(ieee754sp_indef());
661da177e4SLinus Torvalds 		}
671da177e4SLinus Torvalds 		break;
681da177e4SLinus Torvalds 	}
691da177e4SLinus Torvalds 
701da177e4SLinus Torvalds 	ix = x.bits;
711da177e4SLinus Torvalds 
721da177e4SLinus Torvalds 	/* normalize x */
731da177e4SLinus Torvalds 	m = (ix >> 23);
741da177e4SLinus Torvalds 	if (m == 0) {		/* subnormal x */
751da177e4SLinus Torvalds 		for (i = 0; (ix & 0x00800000) == 0; i++)
761da177e4SLinus Torvalds 			ix <<= 1;
771da177e4SLinus Torvalds 		m -= i - 1;
781da177e4SLinus Torvalds 	}
791da177e4SLinus Torvalds 	m -= 127;		/* unbias exponent */
801da177e4SLinus Torvalds 	ix = (ix & 0x007fffff) | 0x00800000;
811da177e4SLinus Torvalds 	if (m & 1)		/* odd m, double x to make it even */
821da177e4SLinus Torvalds 		ix += ix;
831da177e4SLinus Torvalds 	m >>= 1;		/* m = [m/2] */
841da177e4SLinus Torvalds 
851da177e4SLinus Torvalds 	/* generate sqrt(x) bit by bit */
861da177e4SLinus Torvalds 	ix += ix;
871da177e4SLinus Torvalds 	q = s = 0;		/* q = sqrt(x) */
881da177e4SLinus Torvalds 	r = 0x01000000;		/* r = moving bit from right to left */
891da177e4SLinus Torvalds 
901da177e4SLinus Torvalds 	while (r != 0) {
911da177e4SLinus Torvalds 		t = s + r;
921da177e4SLinus Torvalds 		if (t <= ix) {
931da177e4SLinus Torvalds 			s = t + r;
941da177e4SLinus Torvalds 			ix -= t;
951da177e4SLinus Torvalds 			q += r;
961da177e4SLinus Torvalds 		}
971da177e4SLinus Torvalds 		ix += ix;
981da177e4SLinus Torvalds 		r >>= 1;
991da177e4SLinus Torvalds 	}
1001da177e4SLinus Torvalds 
1011da177e4SLinus Torvalds 	if (ix != 0) {
1029e8bad1fSRalf Baechle 		ieee754_setcx(IEEE754_INEXACT);
1031da177e4SLinus Torvalds 		switch (ieee754_csr.rm) {
1041da177e4SLinus Torvalds 		case IEEE754_RP:
1051da177e4SLinus Torvalds 			q += 2;
1061da177e4SLinus Torvalds 			break;
1071da177e4SLinus Torvalds 		case IEEE754_RN:
1081da177e4SLinus Torvalds 			q += (q & 1);
1091da177e4SLinus Torvalds 			break;
1101da177e4SLinus Torvalds 		}
1111da177e4SLinus Torvalds 	}
1121da177e4SLinus Torvalds 	ix = (q >> 1) + 0x3f000000;
1131da177e4SLinus Torvalds 	ix += (m << 23);
1141da177e4SLinus Torvalds 	x.bits = ix;
1151da177e4SLinus Torvalds 	return x;
1161da177e4SLinus Torvalds }
117