1 /* $OpenBSD: dfsqrt.c,v 1.8 2023/03/08 04:43:07 guenther Exp $ */
2 /*
3 (c) Copyright 1986 HEWLETT-PACKARD COMPANY
4 To anyone who acknowledges that this file is provided "AS IS"
5 without any express or implied warranty:
6 permission to use, copy, modify, and distribute this file
7 for any purpose is hereby granted without fee, provided that
8 the above copyright notice and this notice appears in all
9 copies, and that the name of Hewlett-Packard Company not be
10 used in advertising or publicity pertaining to distribution
11 of the software without specific, written prior permission.
12 Hewlett-Packard Company makes no representations about the
13 suitability of this software for any purpose.
14 */
15 /* @(#)dfsqrt.c: Revision: 1.9.88.1 Date: 93/12/07 15:05:46 */
16
17 #include "float.h"
18 #include "dbl_float.h"
19
20 /*
21 * Double Floating-point Square Root
22 */
23
24 int
dbl_fsqrt(srcptr,null,dstptr,status)25 dbl_fsqrt(srcptr, null, dstptr, status)
26 dbl_floating_point *srcptr, *null, *dstptr;
27 unsigned int *status;
28 {
29 register unsigned int srcp1, srcp2, resultp1, resultp2;
30 register unsigned int newbitp1, newbitp2, sump1, sump2;
31 register int src_exponent;
32 register int guardbit = FALSE, even_exponent;
33
34 Dbl_copyfromptr(srcptr,srcp1,srcp2);
35 /*
36 * check source operand for NaN or infinity
37 */
38 if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) {
39 /*
40 * is signaling NaN?
41 */
42 if (Dbl_isone_signaling(srcp1)) {
43 /* trap if INVALIDTRAP enabled */
44 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
45 /* make NaN quiet */
46 Set_invalidflag();
47 Dbl_set_quiet(srcp1);
48 }
49 /*
50 * Return quiet NaN or positive infinity.
51 * Fall thru to negative test if negative infinity.
52 */
53 if (Dbl_iszero_sign(srcp1) ||
54 Dbl_isnotzero_mantissa(srcp1,srcp2)) {
55 Dbl_copytoptr(srcp1,srcp2,dstptr);
56 return(NOEXCEPTION);
57 }
58 }
59
60 /*
61 * check for zero source operand
62 */
63 if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) {
64 Dbl_copytoptr(srcp1,srcp2,dstptr);
65 return(NOEXCEPTION);
66 }
67
68 /*
69 * check for negative source operand
70 */
71 if (Dbl_isone_sign(srcp1)) {
72 /* trap if INVALIDTRAP enabled */
73 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
74 /* make NaN quiet */
75 Set_invalidflag();
76 Dbl_makequietnan(srcp1,srcp2);
77 Dbl_copytoptr(srcp1,srcp2,dstptr);
78 return(NOEXCEPTION);
79 }
80
81 /*
82 * Generate result
83 */
84 if (src_exponent > 0) {
85 even_exponent = Dbl_hidden(srcp1);
86 Dbl_clear_signexponent_set_hidden(srcp1);
87 }
88 else {
89 /* normalize operand */
90 Dbl_clear_signexponent(srcp1);
91 src_exponent++;
92 Dbl_normalize(srcp1,srcp2,src_exponent);
93 even_exponent = src_exponent & 1;
94 }
95 if (even_exponent) {
96 /* exponent is even */
97 /* Add comment here. Explain why odd exponent needs correction */
98 Dbl_leftshiftby1(srcp1,srcp2);
99 }
100 /*
101 * Add comment here. Explain following algorithm.
102 *
103 * Trust me, it works.
104 *
105 */
106 Dbl_setzero(resultp1,resultp2);
107 Dbl_allp1(newbitp1) = 1 << (DBL_P - 32);
108 Dbl_setzero_mantissap2(newbitp2);
109 while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) {
110 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2);
111 if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) {
112 Dbl_leftshiftby1(newbitp1,newbitp2);
113 /* update result */
114 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,
115 resultp1,resultp2);
116 Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2);
117 Dbl_rightshiftby2(newbitp1,newbitp2);
118 }
119 else {
120 Dbl_rightshiftby1(newbitp1,newbitp2);
121 }
122 Dbl_leftshiftby1(srcp1,srcp2);
123 }
124 /* correct exponent for pre-shift */
125 if (even_exponent) {
126 Dbl_rightshiftby1(resultp1,resultp2);
127 }
128
129 /* check for inexact */
130 if (Dbl_isnotzero(srcp1,srcp2)) {
131 if (!even_exponent & Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) {
132 Dbl_increment(resultp1,resultp2);
133 }
134 guardbit = Dbl_lowmantissap2(resultp2);
135 Dbl_rightshiftby1(resultp1,resultp2);
136
137 /* now round result */
138 switch (Rounding_mode()) {
139 case ROUNDPLUS:
140 Dbl_increment(resultp1,resultp2);
141 break;
142 case ROUNDNEAREST:
143 /* stickybit is always true, so guardbit
144 * is enough to determine rounding */
145 if (guardbit) {
146 Dbl_increment(resultp1,resultp2);
147 }
148 break;
149 }
150 /* increment result exponent by 1 if mantissa overflowed */
151 if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2;
152
153 if (Is_inexacttrap_enabled()) {
154 Dbl_set_exponent(resultp1,
155 ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
156 Dbl_copytoptr(resultp1,resultp2,dstptr);
157 return(INEXACTEXCEPTION);
158 }
159 else Set_inexactflag();
160 }
161 else {
162 Dbl_rightshiftby1(resultp1,resultp2);
163 }
164 Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
165 Dbl_copytoptr(resultp1,resultp2,dstptr);
166 return(NOEXCEPTION);
167 }
168