arch/parisc/math-emu/sfsqrt.c

   1 // SPDX-License-Identifier: GPL-2.0-or-later
   2 /*
   3  * Linux/PA-RISC Project (http://www.parisc-linux.org/)
   4  *
   5  * Floating-point emulation code
   6  *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <[email protected]>
   7  */
   8 /*
   9  * BEGIN_DESC
  10  *
  11  *  File:
  12  *      @(#)    pa/spmath/sfsqrt.c              $Revision: 1.1 $
  13  *
  14  *  Purpose:
  15  *      Single Floating-point Square Root
  16  *
  17  *  External Interfaces:
  18  *      sgl_fsqrt(srcptr,nullptr,dstptr,status)
  19  *
  20  *  Internal Interfaces:
  21  *
  22  *  Theory:
  23  *      <<please update with a overview of the operation of this file>>
  24  *
  25  * END_DESC
  26 */
  27
  28
  29 #include "float.h"
  30 #include "sgl_float.h"
  31
  32 /*
  33  *  Single Floating-point Square Root
  34  */
  35
  36 /*ARGSUSED*/
  37 unsigned int
  38 sgl_fsqrt(
  39     sgl_floating_point *srcptr,
  40     unsigned int *nullptr,
  41     sgl_floating_point *dstptr,
  42     unsigned int *status)
  43 {
  44         register unsigned int src, result;
  45         register int src_exponent;
  46         register unsigned int newbit, sum;
  47         register boolean guardbit = FALSE, even_exponent;
  48
  49         src = *srcptr;
  50         /*
  51          * check source operand for NaN or infinity
  52          */
  53         if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) {
  54                 /*
  55                  * is signaling NaN?
  56                  */
  57                 if (Sgl_isone_signaling(src)) {
  58                         /* trap if INVALIDTRAP enabled */
  59                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  60                         /* make NaN quiet */
  61                         Set_invalidflag();
  62                         Sgl_set_quiet(src);
  63                 }
  64                 /*
  65                  * Return quiet NaN or positive infinity.
  66                  *  Fall through to negative test if negative infinity.
  67                  */
  68                 if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) {
  69                         *dstptr = src;
  70                         return(NOEXCEPTION);
  71                 }
  72         }
  73
  74         /*
  75          * check for zero source operand
  76          */
  77         if (Sgl_iszero_exponentmantissa(src)) {
  78                 *dstptr = src;
  79                 return(NOEXCEPTION);
  80         }
  81
  82         /*
  83          * check for negative source operand
  84          */
  85         if (Sgl_isone_sign(src)) {
  86                 /* trap if INVALIDTRAP enabled */
  87                 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  88                 /* make NaN quiet */
  89                 Set_invalidflag();
  90                 Sgl_makequietnan(src);
  91                 *dstptr = src;
  92                 return(NOEXCEPTION);
  93         }
  94
  95         /*
  96          * Generate result
  97          */
  98         if (src_exponent > 0) {
  99                 even_exponent = Sgl_hidden(src);
 100                 Sgl_clear_signexponent_set_hidden(src);
 101         }
 102         else {
 103                 /* normalize operand */
 104                 Sgl_clear_signexponent(src);
 105                 src_exponent++;
 106                 Sgl_normalize(src,src_exponent);
 107                 even_exponent = src_exponent & 1;
 108         }
 109         if (even_exponent) {
 110                 /* exponent is even */
 111                 /* Add comment here.  Explain why odd exponent needs correction */
 112                 Sgl_leftshiftby1(src);
 113         }
 114         /*
 115          * Add comment here.  Explain following algorithm.
 116          *
 117          * Trust me, it works.
 118          *
 119          */
 120         Sgl_setzero(result);
 121         newbit = 1 << SGL_P;
 122         while (newbit && Sgl_isnotzero(src)) {
 123                 Sgl_addition(result,newbit,sum);
 124                 if(sum <= Sgl_all(src)) {
 125                         /* update result */
 126                         Sgl_addition(result,(newbit<<1),result);
 127                         Sgl_subtract(src,sum,src);
 128                 }
 129                 Sgl_rightshiftby1(newbit);
 130                 Sgl_leftshiftby1(src);
 131         }
 132         /* correct exponent for pre-shift */
 133         if (even_exponent) {
 134                 Sgl_rightshiftby1(result);
 135         }
 136
 137         /* check for inexact */
 138         if (Sgl_isnotzero(src)) {
 139                 if (!even_exponent && Sgl_islessthan(result,src))
 140                         Sgl_increment(result);
 141                 guardbit = Sgl_lowmantissa(result);
 142                 Sgl_rightshiftby1(result);
 143
 144                 /*  now round result  */
 145                 switch (Rounding_mode()) {
 146                 case ROUNDPLUS:
 147                      Sgl_increment(result);
 148                      break;
 149                 case ROUNDNEAREST:
 150                      /* stickybit is always true, so guardbit
 151                       * is enough to determine rounding */
 152                      if (guardbit) {
 153                         Sgl_increment(result);
 154                      }
 155                      break;
 156                 }
 157                 /* increment result exponent by 1 if mantissa overflowed */
 158                 if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2;
 159
 160                 if (Is_inexacttrap_enabled()) {
 161                         Sgl_set_exponent(result,
 162                          ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
 163                         *dstptr = result;
 164                         return(INEXACTEXCEPTION);
 165                 }
 166                 else Set_inexactflag();
 167         }
 168         else {
 169                 Sgl_rightshiftby1(result);
 170         }
 171         Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
 172         *dstptr = result;
 173         return(NOEXCEPTION);
 174 }