mips/math-emu/sp_sqrt.c

0001 // SPDX-License-Identifier: GPL-2.0-only
0002 /* IEEE754 floating point arithmetic
0003  * single precision square root
0004  */
0005 /*
0006  * MIPS floating point support
0007  * Copyright (C) 1994-2000 Algorithmics Ltd.
0008  */
0009
0010 #include "ieee754sp.h"
0011
0012 union ieee754sp ieee754sp_sqrt(union ieee754sp x)
0013 {
0014     int ix, s, q, m, t, i;
0015     unsigned int r;
0016     COMPXSP;
0017
0018     /* take care of Inf and NaN */
0019
0020     EXPLODEXSP;
0021     ieee754_clearcx();
0022     FLUSHXSP;
0023
0024     /* x == INF or NAN? */
0025     switch (xc) {
0026     case IEEE754_CLASS_SNAN:
0027         return ieee754sp_nanxcpt(x);
0028
0029     case IEEE754_CLASS_QNAN:
0030         /* sqrt(Nan) = Nan */
0031         return x;
0032
0033     case IEEE754_CLASS_ZERO:
0034         /* sqrt(0) = 0 */
0035         return x;
0036
0037     case IEEE754_CLASS_INF:
0038         if (xs) {
0039             /* sqrt(-Inf) = Nan */
0040             ieee754_setcx(IEEE754_INVALID_OPERATION);
0041             return ieee754sp_indef();
0042         }
0043         /* sqrt(+Inf) = Inf */
0044         return x;
0045
0046     case IEEE754_CLASS_DNORM:
0047     case IEEE754_CLASS_NORM:
0048         if (xs) {
0049             /* sqrt(-x) = Nan */
0050             ieee754_setcx(IEEE754_INVALID_OPERATION);
0051             return ieee754sp_indef();
0052         }
0053         break;
0054     }
0055
0056     ix = x.bits;
0057
0058     /* normalize x */
0059     m = (ix >> 23);
0060     if (m == 0) {       /* subnormal x */
0061         for (i = 0; (ix & 0x00800000) == 0; i++)
0062             ix <<= 1;
0063         m -= i - 1;
0064     }
0065     m -= 127;       /* unbias exponent */
0066     ix = (ix & 0x007fffff) | 0x00800000;
0067     if (m & 1)      /* odd m, double x to make it even */
0068         ix += ix;
0069     m >>= 1;        /* m = [m/2] */
0070
0071     /* generate sqrt(x) bit by bit */
0072     ix += ix;
0073     s = 0;
0074     q = 0;          /* q = sqrt(x) */
0075     r = 0x01000000;     /* r = moving bit from right to left */
0076
0077     while (r != 0) {
0078         t = s + r;
0079         if (t <= ix) {
0080             s = t + r;
0081             ix -= t;
0082             q += r;
0083         }
0084         ix += ix;
0085         r >>= 1;
0086     }
0087
0088     if (ix != 0) {
0089         ieee754_setcx(IEEE754_INEXACT);
0090         switch (ieee754_csr.rm) {
0091         case FPU_CSR_RU:
0092             q += 2;
0093             break;
0094         case FPU_CSR_RN:
0095             q += (q & 1);
0096             break;
0097         }
0098     }
0099     ix = (q >> 1) + 0x3f000000;
0100     ix += (m << 23);
0101     x.bits = ix;
0102     return x;
0103 }