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 |
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 | |