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