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.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 through 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}
196