1 // SPDX-License-Identifier: GPL-2.0-or-later 1
2 /*
3 * Linux/PA-RISC Project (http://www.parisc-li
4 *
5 * Floating-point emulation code
6 * Copyright (C) 2001 Hewlett-Packard (Paul B
7 */
8 /*
9 * BEGIN_DESC
10 *
11 * File:
12 * @(#) pa/spmath/dfsqrt.c
13 *
14 * Purpose:
15 * Double Floating-point Square Root
16 *
17 * External Interfaces:
18 * dbl_fsqrt(srcptr,_nullptr,dstptr,statu
19 *
20 * Internal Interfaces:
21 *
22 * Theory:
23 * <<please update with a overview of the
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
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, re
45 register unsigned int newbitp1, newbit
46 register int src_exponent;
47 register boolean guardbit = FALSE, eve
48
49 Dbl_copyfromptr(srcptr,srcp1,srcp2);
50 /*
51 * check source operand for NaN or inf
52 */
53 if ((src_exponent = Dbl_exponent(srcp1
54 /*
55 * is signaling NaN?
56 */
57 if (Dbl_isone_signaling(srcp1)
58 /* trap if INVALIDTRAP
59 if (Is_invalidtrap_ena
60 /* make NaN quiet */
61 Set_invalidflag();
62 Dbl_set_quiet(srcp1);
63 }
64 /*
65 * Return quiet NaN or positiv
66 * Fall through to negative t
67 */
68 if (Dbl_iszero_sign(srcp1) ||
69 Dbl_isnotzero_mantissa(src
70 Dbl_copytoptr(srcp1,sr
71 return(NOEXCEPTION);
72 }
73 }
74
75 /*
76 * check for zero source operand
77 */
78 if (Dbl_iszero_exponentmantissa(srcp1,
79 Dbl_copytoptr(srcp1,srcp2,dstp
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())
89 /* make NaN quiet */
90 Set_invalidflag();
91 Dbl_makequietnan(srcp1,srcp2);
92 Dbl_copytoptr(srcp1,srcp2,dstp
93 return(NOEXCEPTION);
94 }
95
96 /*
97 * Generate result
98 */
99 if (src_exponent > 0) {
100 even_exponent = Dbl_hidden(src
101 Dbl_clear_signexponent_set_hid
102 }
103 else {
104 /* normalize operand */
105 Dbl_clear_signexponent(srcp1);
106 src_exponent++;
107 Dbl_normalize(srcp1,srcp2,src_
108 even_exponent = src_exponent &
109 }
110 if (even_exponent) {
111 /* exponent is even */
112 /* Add comment here. Explain
113 Dbl_leftshiftby1(srcp1,srcp2);
114 }
115 /*
116 * Add comment here. Explain followin
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
125 Dbl_addition(resultp1,resultp2
126 if(Dbl_isnotgreaterthan(sump1,
127 Dbl_leftshiftby1(newbi
128 /* update result */
129 Dbl_addition(resultp1,
130 resultp1,resultp2);
131 Dbl_subtract(srcp1,src
132 Dbl_rightshiftby2(newb
133 }
134 else {
135 Dbl_rightshiftby1(newb
136 }
137 Dbl_leftshiftby1(srcp1,srcp2);
138 }
139 /* correct exponent for pre-shift */
140 if (even_exponent) {
141 Dbl_rightshiftby1(resultp1,res
142 }
143
144 /* check for inexact */
145 if (Dbl_isnotzero(srcp1,srcp2)) {
146 if (!even_exponent && Dbl_isle
147 Dbl_increment(resultp1
148 }
149 guardbit = Dbl_lowmantissap2(r
150 Dbl_rightshiftby1(resultp1,res
151
152 /* now round result */
153 switch (Rounding_mode()) {
154 case ROUNDPLUS:
155 Dbl_increment(resultp1,re
156 break;
157 case ROUNDNEAREST:
158 /* stickybit is always tr
159 * is enough to determine
160 if (guardbit) {
161 Dbl_increment(resu
162 }
163 break;
164 }
165 /* increment result exponent b
166 if (Dbl_isone_hiddenoverflow(r
167
168 if (Is_inexacttrap_enabled())
169 Dbl_set_exponent(resul
170 ((src_exponent-DBL_BI
171 Dbl_copytoptr(resultp1
172 return(INEXACTEXCEPTIO
173 }
174 else Set_inexactflag();
175 }
176 else {
177 Dbl_rightshiftby1(resultp1,res
178 }
179 Dbl_set_exponent(resultp1,((src_expone
180 Dbl_copytoptr(resultp1,resultp2,dstptr
181 return(NOEXCEPTION);
182 }
183
Linux® is a registered trademark of Linus Torvalds in the United States and other countries.
TOMOYO® is a registered trademark of NTT DATA CORPORATION.