LXR linux/arch/mips/math-emu/sp

   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