sys/arch/hppa/spmath/sfsqrt.c

   1 /*      $NetBSD: sfsqrt.c,v 1.3 2005/12/11 12:17:40 christos Exp $      */
   2
   3 /*      $OpenBSD: sfsqrt.c,v 1.5 2001/03/29 03:58:19 mickey Exp $       */
   4
   5 /*
   6  * Copyright 1996 1995 by Open Software Foundation, Inc.
   7  *              All Rights Reserved
   8  *
   9  * Permission to use, copy, modify, and distribute this software and
  10  * its documentation for any purpose and without fee is hereby granted,
  11  * provided that the above copyright notice appears in all copies and
  12  * that both the copyright notice and this permission notice appear in
  13  * supporting documentation.
  14  *
  15  * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE
  16  * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  17  * FOR A PARTICULAR PURPOSE.
  18  *
  19  * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR
  20  * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
  21  * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT,
  22  * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
  23  * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  24  *
  25  */
  26 /*
  27  * pmk1.1
  28  */
  29 /*
  30  * (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  31  *
  32  * To anyone who acknowledges that this file is provided "AS IS"
  33  * without any express or implied warranty:
  34  *     permission to use, copy, modify, and distribute this file
  35  * for any purpose is hereby granted without fee, provided that
  36  * the above copyright notice and this notice appears in all
  37  * copies, and that the name of Hewlett-Packard Company not be
  38  * used in advertising or publicity pertaining to distribution
  39  * of the software without specific, written prior permission.
  40  * Hewlett-Packard Company makes no representations about the
  41  * suitability of this software for any purpose.
  42  */
  43
  44 #include <sys/cdefs.h>
  45 __KERNEL_RCSID(0, "$NetBSD: sfsqrt.c,v 1.3 2005/12/11 12:17:40 christos Exp $");
  46
  47 #include "../spmath/float.h"
  48 #include "../spmath/sgl_float.h"
  49
  50 /*
  51  *  Single Floating-point Square Root
  52  */
  53
  54 /*ARGSUSED*/
  55 int
  56 sgl_fsqrt(srcptr,dstptr,status)
  57
  58 sgl_floating_point *srcptr, *dstptr;
  59 unsigned int *status;
  60 {
  61         register unsigned int src, result;
  62         register int src_exponent, newbit, sum;
  63         register int guardbit = false, even_exponent;
  64
  65         src = *srcptr;
  66         /*
  67          * check source operand for NaN or infinity
  68          */
  69         if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) {
  70                 /*
  71                  * is signaling NaN?
  72                  */
  73                 if (Sgl_isone_signaling(src)) {
  74                         /* trap if INVALIDTRAP enabled */
  75                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  76                         /* make NaN quiet */
  77                         Set_invalidflag();
  78                         Sgl_set_quiet(src);
  79                 }
  80                 /*
  81                  * Return quiet NaN or positive infinity.
  82                  *  Fall thru to negative test if negative infinity.
  83                  */
  84                 if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) {
  85                         *dstptr = src;
  86                         return(NOEXCEPTION);
  87                 }
  88         }
  89
  90         /*
  91          * check for zero source operand
  92          */
  93         if (Sgl_iszero_exponentmantissa(src)) {
  94                 *dstptr = src;
  95                 return(NOEXCEPTION);
  96         }
  97
  98         /*
  99          * check for negative source operand
 100          */
 101         if (Sgl_isone_sign(src)) {
 102                 /* trap if INVALIDTRAP enabled */
 103                 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 104                 /* make NaN quiet */
 105                 Set_invalidflag();
 106                 Sgl_makequietnan(src);
 107                 *dstptr = src;
 108                 return(NOEXCEPTION);
 109         }
 110
 111         /*
 112          * Generate result
 113          */
 114         if (src_exponent > 0) {
 115                 even_exponent = Sgl_hidden(src);
 116                 Sgl_clear_signexponent_set_hidden(src);
 117         }
 118         else {
 119                 /* normalize operand */
 120                 Sgl_clear_signexponent(src);
 121                 src_exponent++;
 122                 Sgl_normalize(src,src_exponent);
 123                 even_exponent = src_exponent & 1;
 124         }
 125         if (even_exponent) {
 126                 /* exponent is even */
 127                 /* Add comment here.  Explain why odd exponent needs correction */
 128                 Sgl_leftshiftby1(src);
 129         }
 130         /*
 131          * Add comment here.  Explain following algorithm.
 132          *
 133          * Trust me, it works.
 134          *
 135          */
 136         Sgl_setzero(result);
 137         newbit = 1 << SGL_P;
 138         while (newbit && Sgl_isnotzero(src)) {
 139                 Sgl_addition(result,newbit,sum);
 140                 if(sum <= Sgl_all(src)) {
 141                         /* update result */
 142                         Sgl_addition(result,(newbit<<1),result);
 143                         Sgl_subtract(src,sum,src);
 144                 }
 145                 Sgl_rightshiftby1(newbit);
 146                 Sgl_leftshiftby1(src);
 147         }
 148         /* correct exponent for pre-shift */
 149         if (even_exponent) {
 150                 Sgl_rightshiftby1(result);
 151         }
 152
 153         /* check for inexact */
 154         if (Sgl_isnotzero(src)) {
 155                 if (!even_exponent & Sgl_islessthan(result,src))
 156                         Sgl_increment(result);
 157                 guardbit = Sgl_lowmantissa(result);
 158                 Sgl_rightshiftby1(result);
 159
 160                 /*  now round result  */
 161                 switch (Rounding_mode()) {
 162                 case ROUNDPLUS:
 163                      Sgl_increment(result);
 164                      break;
 165                 case ROUNDNEAREST:
 166                      /* stickybit is always true, so guardbit
 167                       * is enough to determine rounding */
 168                      if (guardbit) {
 169                         Sgl_increment(result);
 170                      }
 171                      break;
 172                 }
 173                 /* increment result exponent by 1 if mantissa overflowed */
 174                 if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2;
 175
 176                 if (Is_inexacttrap_enabled()) {
 177                         Sgl_set_exponent(result,
 178                          ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
 179                         *dstptr = result;
 180                         return(INEXACTEXCEPTION);
 181                 }
 182                 else Set_inexactflag();
 183         }
 184         else {
 185                 Sgl_rightshiftby1(result);
 186         }
 187         Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
 188         *dstptr = result;
 189         return(NOEXCEPTION);
 190 }