Quelle  sp_sqrt.c   Sprache: C

// SPDX-License-Identifier: GPL-2.0-only
/* IEEE754 floating point arithmetic
 * single precision square root
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.
 */

#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;
 s = 0;
 q = 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;
}