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/sfsqrt.c $Revision: 1.1 $ |
13 | * |
14 | * Purpose: |
15 | * Single Floating-point Square Root |
16 | * |
17 | * External Interfaces: |
18 | * sgl_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 "sgl_float.h" |
31 | |
32 | /* |
33 | * Single Floating-point Square Root |
34 | */ |
35 | |
36 | /*ARGSUSED*/ |
37 | unsigned int |
38 | sgl_fsqrt( |
39 | sgl_floating_point *srcptr, |
40 | unsigned int *_nullptr, |
41 | sgl_floating_point *dstptr, |
42 | unsigned int *status) |
43 | { |
44 | register unsigned int src, result; |
45 | register int src_exponent; |
46 | register unsigned int newbit, sum; |
47 | register boolean guardbit = FALSE, even_exponent; |
48 | |
49 | src = *srcptr; |
50 | /* |
51 | * check source operand for NaN or infinity |
52 | */ |
53 | if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) { |
54 | /* |
55 | * is signaling NaN? |
56 | */ |
57 | if (Sgl_isone_signaling(src)) { |
58 | /* trap if INVALIDTRAP enabled */ |
59 | if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); |
60 | /* make NaN quiet */ |
61 | Set_invalidflag(); |
62 | Sgl_set_quiet(src); |
63 | } |
64 | /* |
65 | * Return quiet NaN or positive infinity. |
66 | * Fall through to negative test if negative infinity. |
67 | */ |
68 | if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) { |
69 | *dstptr = src; |
70 | return(NOEXCEPTION); |
71 | } |
72 | } |
73 | |
74 | /* |
75 | * check for zero source operand |
76 | */ |
77 | if (Sgl_iszero_exponentmantissa(src)) { |
78 | *dstptr = src; |
79 | return(NOEXCEPTION); |
80 | } |
81 | |
82 | /* |
83 | * check for negative source operand |
84 | */ |
85 | if (Sgl_isone_sign(src)) { |
86 | /* trap if INVALIDTRAP enabled */ |
87 | if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); |
88 | /* make NaN quiet */ |
89 | Set_invalidflag(); |
90 | Sgl_makequietnan(src); |
91 | *dstptr = src; |
92 | return(NOEXCEPTION); |
93 | } |
94 | |
95 | /* |
96 | * Generate result |
97 | */ |
98 | if (src_exponent > 0) { |
99 | even_exponent = Sgl_hidden(src); |
100 | Sgl_clear_signexponent_set_hidden(src); |
101 | } |
102 | else { |
103 | /* normalize operand */ |
104 | Sgl_clear_signexponent(src); |
105 | src_exponent++; |
106 | Sgl_normalize(src,src_exponent); |
107 | even_exponent = src_exponent & 1; |
108 | } |
109 | if (even_exponent) { |
110 | /* exponent is even */ |
111 | /* Add comment here. Explain why odd exponent needs correction */ |
112 | Sgl_leftshiftby1(src); |
113 | } |
114 | /* |
115 | * Add comment here. Explain following algorithm. |
116 | * |
117 | * Trust me, it works. |
118 | * |
119 | */ |
120 | Sgl_setzero(result); |
121 | newbit = 1 << SGL_P; |
122 | while (newbit && Sgl_isnotzero(src)) { |
123 | Sgl_addition(result,newbit,sum); |
124 | if(sum <= Sgl_all(src)) { |
125 | /* update result */ |
126 | Sgl_addition(result,(newbit<<1),result); |
127 | Sgl_subtract(src,sum,src); |
128 | } |
129 | Sgl_rightshiftby1(newbit); |
130 | Sgl_leftshiftby1(src); |
131 | } |
132 | /* correct exponent for pre-shift */ |
133 | if (even_exponent) { |
134 | Sgl_rightshiftby1(result); |
135 | } |
136 | |
137 | /* check for inexact */ |
138 | if (Sgl_isnotzero(src)) { |
139 | if (!even_exponent && Sgl_islessthan(result,src)) |
140 | Sgl_increment(result); |
141 | guardbit = Sgl_lowmantissa(result); |
142 | Sgl_rightshiftby1(result); |
143 | |
144 | /* now round result */ |
145 | switch (Rounding_mode()) { |
146 | case ROUNDPLUS: |
147 | Sgl_increment(result); |
148 | break; |
149 | case ROUNDNEAREST: |
150 | /* stickybit is always true, so guardbit |
151 | * is enough to determine rounding */ |
152 | if (guardbit) { |
153 | Sgl_increment(result); |
154 | } |
155 | break; |
156 | } |
157 | /* increment result exponent by 1 if mantissa overflowed */ |
158 | if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2; |
159 | |
160 | if (Is_inexacttrap_enabled()) { |
161 | Sgl_set_exponent(result, |
162 | ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); |
163 | *dstptr = result; |
164 | return(INEXACTEXCEPTION); |
165 | } |
166 | else Set_inexactflag(); |
167 | } |
168 | else { |
169 | Sgl_rightshiftby1(result); |
170 | } |
171 | Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); |
172 | *dstptr = result; |
173 | return(NOEXCEPTION); |
174 | } |
175 | |