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