sys/arch/hppa/spmath/dfsqrt.c

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