arch/parisc/math-emu/dfsqrt.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) <bame@debian.org>
   7  */
   8 /*
   9  * BEGIN_DESC
  10  *
  11  *  File:
  12  *      @(#)    pa/spmath/dfsqrt.c              $Revision: 1.1 $
  13  *
  14  *  Purpose:
  15  *      Double Floating-point Square Root
  16  *
  17  *  External Interfaces:
  18  *      dbl_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 "dbl_float.h"
  31
  32 /*
  33  *  Double Floating-point Square Root
  34  */
  35
  36 /*ARGSUSED*/
  37 unsigned int
  38 dbl_fsqrt(
  39             dbl_floating_point *srcptr,
  40             unsigned int *nullptr,
  41             dbl_floating_point *dstptr,
  42             unsigned int *status)
  43 {
  44         register unsigned int srcp1, srcp2, resultp1, resultp2;
  45         register unsigned int newbitp1, newbitp2, sump1, sump2;
  46         register int src_exponent;
  47         register boolean guardbit = FALSE, even_exponent;
  48
  49         Dbl_copyfromptr(srcptr,srcp1,srcp2);
  50         /*
  51          * check source operand for NaN or infinity
  52          */
  53         if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) {
  54                 /*
  55                  * is signaling NaN?
  56                  */
  57                 if (Dbl_isone_signaling(srcp1)) {
  58                         /* trap if INVALIDTRAP enabled */
  59                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  60                         /* make NaN quiet */
  61                         Set_invalidflag();
  62                         Dbl_set_quiet(srcp1);
  63                 }
  64                 /*
  65                  * Return quiet NaN or positive infinity.
  66                  *  Fall through to negative test if negative infinity.
  67                  */
  68                 if (Dbl_iszero_sign(srcp1) ||
  69                     Dbl_isnotzero_mantissa(srcp1,srcp2)) {
  70                         Dbl_copytoptr(srcp1,srcp2,dstptr);
  71                         return(NOEXCEPTION);
  72                 }
  73         }
  74
  75         /*
  76          * check for zero source operand
  77          */
  78         if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) {
  79                 Dbl_copytoptr(srcp1,srcp2,dstptr);
  80                 return(NOEXCEPTION);
  81         }
  82
  83         /*
  84          * check for negative source operand
  85          */
  86         if (Dbl_isone_sign(srcp1)) {
  87                 /* trap if INVALIDTRAP enabled */
  88                 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  89                 /* make NaN quiet */
  90                 Set_invalidflag();
  91                 Dbl_makequietnan(srcp1,srcp2);
  92                 Dbl_copytoptr(srcp1,srcp2,dstptr);
  93                 return(NOEXCEPTION);
  94         }
  95
  96         /*
  97          * Generate result
  98          */
  99         if (src_exponent > 0) {
 100                 even_exponent = Dbl_hidden(srcp1);
 101                 Dbl_clear_signexponent_set_hidden(srcp1);
 102         }
 103         else {
 104                 /* normalize operand */
 105                 Dbl_clear_signexponent(srcp1);
 106                 src_exponent++;
 107                 Dbl_normalize(srcp1,srcp2,src_exponent);
 108                 even_exponent = src_exponent & 1;
 109         }
 110         if (even_exponent) {
 111                 /* exponent is even */
 112                 /* Add comment here.  Explain why odd exponent needs correction */
 113                 Dbl_leftshiftby1(srcp1,srcp2);
 114         }
 115         /*
 116          * Add comment here.  Explain following algorithm.
 117          *
 118          * Trust me, it works.
 119          *
 120          */
 121         Dbl_setzero(resultp1,resultp2);
 122         Dbl_allp1(newbitp1) = 1 << (DBL_P - 32);
 123         Dbl_setzero_mantissap2(newbitp2);
 124         while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) {
 125                 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2);
 126                 if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) {
 127                         Dbl_leftshiftby1(newbitp1,newbitp2);
 128                         /* update result */
 129                         Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,
 130                          resultp1,resultp2);
 131                         Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2);
 132                         Dbl_rightshiftby2(newbitp1,newbitp2);
 133                 }
 134                 else {
 135                         Dbl_rightshiftby1(newbitp1,newbitp2);
 136                 }
 137                 Dbl_leftshiftby1(srcp1,srcp2);
 138         }
 139         /* correct exponent for pre-shift */
 140         if (even_exponent) {
 141                 Dbl_rightshiftby1(resultp1,resultp2);
 142         }
 143
 144         /* check for inexact */
 145         if (Dbl_isnotzero(srcp1,srcp2)) {
 146                 if (!even_exponent && Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) {
 147                         Dbl_increment(resultp1,resultp2);
 148                 }
 149                 guardbit = Dbl_lowmantissap2(resultp2);
 150                 Dbl_rightshiftby1(resultp1,resultp2);
 151
 152                 /*  now round result  */
 153                 switch (Rounding_mode()) {
 154                 case ROUNDPLUS:
 155                      Dbl_increment(resultp1,resultp2);
 156                      break;
 157                 case ROUNDNEAREST:
 158                      /* stickybit is always true, so guardbit
 159                       * is enough to determine rounding */
 160                      if (guardbit) {
 161                             Dbl_increment(resultp1,resultp2);
 162                      }
 163                      break;
 164                 }
 165                 /* increment result exponent by 1 if mantissa overflowed */
 166                 if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2;
 167
 168                 if (Is_inexacttrap_enabled()) {
 169                         Dbl_set_exponent(resultp1,
 170                          ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
 171                         Dbl_copytoptr(resultp1,resultp2,dstptr);
 172                         return(INEXACTEXCEPTION);
 173                 }
 174                 else Set_inexactflag();
 175         }
 176         else {
 177                 Dbl_rightshiftby1(resultp1,resultp2);
 178         }
 179         Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
 180         Dbl_copytoptr(resultp1,resultp2,dstptr);
 181         return(NOEXCEPTION);
 182 }