~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~

TOMOYO Linux Cross Reference
Linux/arch/parisc/math-emu/sfsqrt.c

Version: ~ [ linux-5.6 ] ~ [ linux-5.5.13 ] ~ [ linux-5.4.28 ] ~ [ linux-5.3.18 ] ~ [ linux-5.2.21 ] ~ [ linux-5.1.21 ] ~ [ linux-5.0.21 ] ~ [ linux-4.20.17 ] ~ [ linux-4.19.113 ] ~ [ linux-4.18.20 ] ~ [ linux-4.17.19 ] ~ [ linux-4.16.18 ] ~ [ linux-4.15.18 ] ~ [ linux-4.14.174 ] ~ [ linux-4.13.16 ] ~ [ linux-4.12.14 ] ~ [ linux-4.11.12 ] ~ [ linux-4.10.17 ] ~ [ linux-4.9.217 ] ~ [ linux-4.8.17 ] ~ [ linux-4.7.10 ] ~ [ linux-4.6.7 ] ~ [ linux-4.5.7 ] ~ [ linux-4.4.217 ] ~ [ linux-4.3.6 ] ~ [ linux-4.2.8 ] ~ [ linux-4.1.52 ] ~ [ linux-4.0.9 ] ~ [ linux-3.19.8 ] ~ [ linux-3.18.140 ] ~ [ linux-3.17.8 ] ~ [ linux-3.16.82 ] ~ [ linux-3.15.10 ] ~ [ linux-3.14.79 ] ~ [ linux-3.13.11 ] ~ [ linux-3.12.74 ] ~ [ linux-3.11.10 ] ~ [ linux-3.10.108 ] ~ [ linux-3.9.11 ] ~ [ linux-3.8.13 ] ~ [ linux-3.7.10 ] ~ [ linux-3.6.11 ] ~ [ linux-3.5.7 ] ~ [ linux-3.4.113 ] ~ [ linux-3.3.8 ] ~ [ linux-3.2.102 ] ~ [ linux-3.1.10 ] ~ [ linux-3.0.101 ] ~ [ linux-2.6.32.71 ] ~ [ linux-2.6.0 ] ~ [ linux-2.4.37.11 ] ~ [ unix-v6-master ] ~ [ ccs-tools-1.8.5 ] ~ [ policy-sample ] ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

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

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~

kernel.org | git.kernel.org | LWN.net | Project Home | Wiki (Japanese) | Wiki (English) | SVN repository | Mail admin

Linux® is a registered trademark of Linus Torvalds in the United States and other countries.
TOMOYO® is a registered trademark of NTT DATA CORPORATION.

osdn.jp