arch/mips/math-emu/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
  12 union 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 }