123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114 |
- /* IEEE754 floating point arithmetic
- * single precision square root
- */
- /*
- * MIPS floating point support
- * Copyright (C) 1994-2000 Algorithmics Ltd.
- *
- * This program is free software; you can distribute it and/or modify it
- * under the terms of the GNU General Public License (Version 2) as
- * published by the Free Software Foundation.
- *
- * This program is distributed in the hope it will be useful, but WITHOUT
- * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- * for more details.
- *
- * You should have received a copy of the GNU General Public License along
- * with this program; if not, write to the Free Software Foundation, Inc.,
- * 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
- */
- #include "ieee754sp.h"
- union ieee754sp ieee754sp_sqrt(union ieee754sp x)
- {
- int ix, s, q, m, t, i;
- unsigned int r;
- COMPXSP;
- /* take care of Inf and NaN */
- EXPLODEXSP;
- ieee754_clearcx();
- FLUSHXSP;
- /* x == INF or NAN? */
- switch (xc) {
- case IEEE754_CLASS_SNAN:
- return ieee754sp_nanxcpt(x);
- case IEEE754_CLASS_QNAN:
- /* sqrt(Nan) = Nan */
- return x;
- case IEEE754_CLASS_ZERO:
- /* sqrt(0) = 0 */
- return x;
- case IEEE754_CLASS_INF:
- if (xs) {
- /* sqrt(-Inf) = Nan */
- ieee754_setcx(IEEE754_INVALID_OPERATION);
- return ieee754sp_indef();
- }
- /* sqrt(+Inf) = Inf */
- return x;
- case IEEE754_CLASS_DNORM:
- case IEEE754_CLASS_NORM:
- if (xs) {
- /* sqrt(-x) = Nan */
- ieee754_setcx(IEEE754_INVALID_OPERATION);
- return ieee754sp_indef();
- }
- break;
- }
- ix = x.bits;
- /* normalize x */
- m = (ix >> 23);
- if (m == 0) { /* subnormal x */
- for (i = 0; (ix & 0x00800000) == 0; i++)
- ix <<= 1;
- m -= i - 1;
- }
- m -= 127; /* unbias exponent */
- ix = (ix & 0x007fffff) | 0x00800000;
- if (m & 1) /* odd m, double x to make it even */
- ix += ix;
- m >>= 1; /* m = [m/2] */
- /* generate sqrt(x) bit by bit */
- ix += ix;
- q = s = 0; /* q = sqrt(x) */
- r = 0x01000000; /* r = moving bit from right to left */
- while (r != 0) {
- t = s + r;
- if (t <= ix) {
- s = t + r;
- ix -= t;
- q += r;
- }
- ix += ix;
- r >>= 1;
- }
- if (ix != 0) {
- ieee754_setcx(IEEE754_INEXACT);
- switch (ieee754_csr.rm) {
- case FPU_CSR_RU:
- q += 2;
- break;
- case FPU_CSR_RN:
- q += (q & 1);
- break;
- }
- }
- ix = (q >> 1) + 0x3f000000;
- ix += (m << 23);
- x.bits = ix;
- return x;
- }
|