1 // SPDX-License-Identifier: GPL-2.0-or-later
2 /*
3  * Linux/PA-RISC Project (http://www.parisc-linux.org/)
4  *
5  * Floating-point emulation code
6  *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
7  */
8 /*
9  * BEGIN_DESC
10  *
11  *  File:
12  *	@(#)	pa/spmath/dfsqrt.c		$Revision: 1.1 $
13  *
14  *  Purpose:
15  *	Double Floating-point Square Root
16  *
17  *  External Interfaces:
18  *	dbl_fsqrt(srcptr,nullptr,dstptr,status)
19  *
20  *  Internal Interfaces:
21  *
22  *  Theory:
23  *	<<please update with a overview of the operation of this file>>
24  *
25  * END_DESC
26 */
27 
28 
29 #include "float.h"
30 #include "dbl_float.h"
31 
32 /*
33  *  Double Floating-point Square Root
34  */
35 
36 /*ARGSUSED*/
37 unsigned int
dbl_fsqrt(dbl_floating_point * srcptr,unsigned int * nullptr,dbl_floating_point * dstptr,unsigned int * status)38 dbl_fsqrt(
39 	    dbl_floating_point *srcptr,
40 	    unsigned int *nullptr,
41 	    dbl_floating_point *dstptr,
42 	    unsigned int *status)
43 {
44 	register unsigned int srcp1, srcp2, resultp1, resultp2;
45 	register unsigned int newbitp1, newbitp2, sump1, sump2;
46 	register int src_exponent;
47 	register boolean guardbit = FALSE, even_exponent;
48 
49 	Dbl_copyfromptr(srcptr,srcp1,srcp2);
50         /*
51          * check source operand for NaN or infinity
52          */
53         if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) {
54                 /*
55                  * is signaling NaN?
56                  */
57                 if (Dbl_isone_signaling(srcp1)) {
58                         /* trap if INVALIDTRAP enabled */
59                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
60                         /* make NaN quiet */
61                         Set_invalidflag();
62                         Dbl_set_quiet(srcp1);
63                 }
64                 /*
65                  * Return quiet NaN or positive infinity.
66 		 *  Fall through to negative test if negative infinity.
67                  */
68 		if (Dbl_iszero_sign(srcp1) ||
69 		    Dbl_isnotzero_mantissa(srcp1,srcp2)) {
70                 	Dbl_copytoptr(srcp1,srcp2,dstptr);
71                 	return(NOEXCEPTION);
72 		}
73         }
74 
75         /*
76          * check for zero source operand
77          */
78 	if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) {
79 		Dbl_copytoptr(srcp1,srcp2,dstptr);
80 		return(NOEXCEPTION);
81 	}
82 
83         /*
84          * check for negative source operand
85          */
86 	if (Dbl_isone_sign(srcp1)) {
87 		/* trap if INVALIDTRAP enabled */
88 		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
89 		/* make NaN quiet */
90 		Set_invalidflag();
91 		Dbl_makequietnan(srcp1,srcp2);
92 		Dbl_copytoptr(srcp1,srcp2,dstptr);
93 		return(NOEXCEPTION);
94 	}
95 
96 	/*
97 	 * Generate result
98 	 */
99 	if (src_exponent > 0) {
100 		even_exponent = Dbl_hidden(srcp1);
101 		Dbl_clear_signexponent_set_hidden(srcp1);
102 	}
103 	else {
104 		/* normalize operand */
105 		Dbl_clear_signexponent(srcp1);
106 		src_exponent++;
107 		Dbl_normalize(srcp1,srcp2,src_exponent);
108 		even_exponent = src_exponent & 1;
109 	}
110 	if (even_exponent) {
111 		/* exponent is even */
112 		/* Add comment here.  Explain why odd exponent needs correction */
113 		Dbl_leftshiftby1(srcp1,srcp2);
114 	}
115 	/*
116 	 * Add comment here.  Explain following algorithm.
117 	 *
118 	 * Trust me, it works.
119 	 *
120 	 */
121 	Dbl_setzero(resultp1,resultp2);
122 	Dbl_allp1(newbitp1) = 1 << (DBL_P - 32);
123 	Dbl_setzero_mantissap2(newbitp2);
124 	while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) {
125 		Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2);
126 		if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) {
127 			Dbl_leftshiftby1(newbitp1,newbitp2);
128 			/* update result */
129 			Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,
130 			 resultp1,resultp2);
131 			Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2);
132 			Dbl_rightshiftby2(newbitp1,newbitp2);
133 		}
134 		else {
135 			Dbl_rightshiftby1(newbitp1,newbitp2);
136 		}
137 		Dbl_leftshiftby1(srcp1,srcp2);
138 	}
139 	/* correct exponent for pre-shift */
140 	if (even_exponent) {
141 		Dbl_rightshiftby1(resultp1,resultp2);
142 	}
143 
144 	/* check for inexact */
145 	if (Dbl_isnotzero(srcp1,srcp2)) {
146 		if (!even_exponent && Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) {
147 			Dbl_increment(resultp1,resultp2);
148 		}
149 		guardbit = Dbl_lowmantissap2(resultp2);
150 		Dbl_rightshiftby1(resultp1,resultp2);
151 
152 		/*  now round result  */
153 		switch (Rounding_mode()) {
154 		case ROUNDPLUS:
155 		     Dbl_increment(resultp1,resultp2);
156 		     break;
157 		case ROUNDNEAREST:
158 		     /* stickybit is always true, so guardbit
159 		      * is enough to determine rounding */
160 		     if (guardbit) {
161 			    Dbl_increment(resultp1,resultp2);
162 		     }
163 		     break;
164 		}
165 		/* increment result exponent by 1 if mantissa overflowed */
166 		if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2;
167 
168 		if (Is_inexacttrap_enabled()) {
169 			Dbl_set_exponent(resultp1,
170 			 ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
171 			Dbl_copytoptr(resultp1,resultp2,dstptr);
172 			return(INEXACTEXCEPTION);
173 		}
174 		else Set_inexactflag();
175 	}
176 	else {
177 		Dbl_rightshiftby1(resultp1,resultp2);
178 	}
179 	Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
180 	Dbl_copytoptr(resultp1,resultp2,dstptr);
181 	return(NOEXCEPTION);
182 }
183