PageRenderTime 30ms CodeModel.GetById 19ms app.highlight 8ms RepoModel.GetById 1ms app.codeStats 0ms

/arch/parisc/math-emu/dfsqrt.c

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