Audacity 3.2.0
Namespaces | Functions
simd_complex_conversions Namespace Reference

Namespaces

namespace  details
 

Functions

__m128 atan_ps (__m128 x)
 
__m128 atan2_ps (__m128 y, __m128 x)
 
std::pair< __m128, __m128 > sincos_ps (__m128 x)
 
float atan2_ss (float y, float x)
 
std::pair< float, float > sincos_ss (float angle)
 
__m128 norm (__m128 x, __m128 y)
 
float sqrt_ss (float x)
 
template<typename fnc >
void perform_parallel_simd_aligned (const std::complex< float > *input, float *output, int n, const fnc &f)
 
void rotate_parallel_simd_aligned (const float *oldPhase, const float *newPhase, std::complex< float > *output, int n)
 

Function Documentation

◆ atan2_ps()

__m128 simd_complex_conversions::atan2_ps ( __m128  y,
__m128  x 
)
inline

Definition at line 129 of file SimdComplexConversions_sse2.h.

130{
131 using namespace details;
132
133 __m128 zero = _mm_setzero_ps();
134 __m128 x_eq_0 = _mm_cmpeq_ps(x, zero);
135 __m128 x_gt_0 = _mm_cmpgt_ps(x, zero);
136 __m128 x_le_0 = _mm_cmple_ps(x, zero);
137 __m128 y_eq_0 = _mm_cmpeq_ps(y, zero);
138 __m128 x_lt_0 = _mm_cmplt_ps(x, zero);
139 __m128 y_lt_0 = _mm_cmplt_ps(y, zero);
140
141 __m128 zero_mask = _mm_and_ps(x_eq_0, y_eq_0);
142 __m128 zero_mask_other_case = _mm_and_ps(y_eq_0, x_gt_0);
143 zero_mask = _mm_or_ps(zero_mask, zero_mask_other_case);
144
145 __m128 pio2_mask = _mm_andnot_ps(y_eq_0, x_eq_0);
146 __m128 pio2_mask_sign = _mm_and_ps(y_lt_0, _mm_set1_ps(sign_mask));
147 __m128 pio2_result = _mm_set1_ps(cephes_PIO2F);
148 pio2_result = _mm_xor_ps(pio2_result, pio2_mask_sign);
149 pio2_result = _mm_and_ps(pio2_mask, pio2_result);
150
151 __m128 pi_mask = _mm_and_ps(y_eq_0, x_lt_0);
152 __m128 pi = _mm_set1_ps(cephes_PIF);
153 __m128 pi_result = _mm_and_ps(pi_mask, pi);
154
155 __m128 swap_sign_mask_offset = _mm_and_ps(x_lt_0, y_lt_0);
156 swap_sign_mask_offset =
157 _mm_and_ps(swap_sign_mask_offset, _mm_set1_ps(sign_mask));
158
159 __m128 offset0 = _mm_setzero_ps();
160 __m128 offset1 = _mm_set1_ps(cephes_PIF);
161 offset1 = _mm_xor_ps(offset1, swap_sign_mask_offset);
162
163 __m128 offset = _mm_andnot_ps(x_lt_0, offset0);
164 offset = _mm_and_ps(x_lt_0, offset1);
165
166 __m128 arg = _mm_div_ps(y, x);
167 __m128 atan_result = atan_ps(arg);
168 atan_result = _mm_add_ps(atan_result, offset);
169
170 /* select between zero_result, pio2_result and atan_result */
171
172 __m128 result = _mm_andnot_ps(zero_mask, pio2_result);
173 atan_result = _mm_andnot_ps(zero_mask, atan_result);
174 atan_result = _mm_andnot_ps(pio2_mask, atan_result);
175 result = _mm_or_ps(result, atan_result);
176 result = _mm_or_ps(result, pi_result);
177
178 return result;
179}

References atan_ps(), simd_complex_conversions::details::cephes_PIF, simd_complex_conversions::details::cephes_PIO2F, MIR::anonymous_namespace{MirUtils.cpp}::pi, and simd_complex_conversions::details::sign_mask.

Referenced by atan2_ss().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ atan2_ss()

float simd_complex_conversions::atan2_ss ( float  y,
float  x 
)
inline

Definition at line 275 of file SimdComplexConversions_sse2.h.

276{
277 return _mm_cvtss_f32(atan2_ps(_mm_set_ss(y), _mm_set_ss(x)));
278}
__m128 atan2_ps(__m128 y, __m128 x)

References atan2_ps().

Here is the call graph for this function:

◆ atan_ps()

__m128 simd_complex_conversions::atan_ps ( __m128  x)
inline

Definition at line 68 of file SimdComplexConversions_sse2.h.

69{
70 using namespace details;
71
72 __m128 sign_bit, y;
73
74 sign_bit = x;
75 /* take the absolute value */
76 x = _mm_and_ps(x, _mm_set1_ps(inv_sign_mask));
77 /* extract the sign bit (upper one) */
78 sign_bit = _mm_and_ps(sign_bit, _mm_set1_ps(sign_mask));
79
80 /* range reduction, init x and y depending on range */
81 /* x > 2.414213562373095 */
82 __m128 cmp0 = _mm_cmpgt_ps(x, _mm_set1_ps(2.414213562373095f));
83 /* x > 0.4142135623730950 */
84 __m128 cmp1 = _mm_cmpgt_ps(x, _mm_set1_ps(0.4142135623730950f));
85
86 /* x > 0.4142135623730950 && !( x > 2.414213562373095 ) */
87 __m128 cmp2 = _mm_andnot_ps(cmp0, cmp1);
88
89 /* -( 1.0/x ) */
90 __m128 y0 = _mm_and_ps(cmp0, _mm_set1_ps(cephes_PIO2F));
91 __m128 x0 = _mm_div_ps(_mm_set1_ps(1.0f), x);
92 x0 = _mm_xor_ps(x0, _mm_set1_ps(sign_mask));
93
94 __m128 y1 = _mm_and_ps(cmp2, _mm_set1_ps(cephes_PIO4F));
95 /* (x-1.0)/(x+1.0) */
96 __m128 x1_o = _mm_sub_ps(x, _mm_set1_ps(1.0f));
97 __m128 x1_u = _mm_add_ps(x, _mm_set1_ps(1.0f));
98 __m128 x1 = _mm_div_ps(x1_o, x1_u);
99
100 __m128 x2 = _mm_and_ps(cmp2, x1);
101 x0 = _mm_and_ps(cmp0, x0);
102 x2 = _mm_or_ps(x2, x0);
103 cmp1 = _mm_or_ps(cmp0, cmp2);
104 x2 = _mm_and_ps(cmp1, x2);
105 x = _mm_andnot_ps(cmp1, x);
106 x = _mm_or_ps(x2, x);
107
108 y = _mm_or_ps(y0, y1);
109
110 __m128 zz = _mm_mul_ps(x, x);
111 __m128 acc = _mm_set1_ps(atancof_p0);
112 acc = _mm_mul_ps(acc, zz);
113 acc = _mm_sub_ps(acc, _mm_set1_ps(atancof_p1));
114 acc = _mm_mul_ps(acc, zz);
115 acc = _mm_add_ps(acc, _mm_set1_ps(atancof_p2));
116 acc = _mm_mul_ps(acc, zz);
117 acc = _mm_sub_ps(acc, _mm_set1_ps(atancof_p3));
118 acc = _mm_mul_ps(acc, zz);
119 acc = _mm_mul_ps(acc, x);
120 acc = _mm_add_ps(acc, x);
121 y = _mm_add_ps(y, acc);
122
123 /* update the sign */
124 y = _mm_xor_ps(y, sign_bit);
125
126 return y;
127}

References simd_complex_conversions::details::atancof_p0, simd_complex_conversions::details::atancof_p1, simd_complex_conversions::details::atancof_p2, simd_complex_conversions::details::atancof_p3, simd_complex_conversions::details::cephes_PIO2F, simd_complex_conversions::details::cephes_PIO4F, simd_complex_conversions::details::inv_sign_mask, and simd_complex_conversions::details::sign_mask.

Referenced by atan2_ps().

Here is the caller graph for this function:

◆ norm()

__m128 simd_complex_conversions::norm ( __m128  x,
__m128  y 
)
inline

Definition at line 286 of file SimdComplexConversions_sse2.h.

287{
288 return _mm_add_ps(_mm_mul_ps(x, x), _mm_mul_ps(y, y));
289}

Referenced by staffpad::vo::calcNorms(), and graphics::Normalized().

Here is the caller graph for this function:

◆ perform_parallel_simd_aligned()

template<typename fnc >
void simd_complex_conversions::perform_parallel_simd_aligned ( const std::complex< float > *  input,
float *  output,
int  n,
const fnc &  f 
)

Definition at line 299 of file SimdComplexConversions_sse2.h.

301{
302 for (int i = 0; i <= n - 4; i += 4)
303 {
304 // Safe according to C++ standard
305 auto p1 = _mm_load_ps(reinterpret_cast<const float*>(input + i));
306 auto p2 = _mm_load_ps(reinterpret_cast<const float*>(input + i + 2));
307
308 // p1 = {real(c1), imag(c1), real(c2), imag(c2)}
309 // p2 = {real(c3), imag(c3), real(c4), imag(c4)}
310
311 auto rp = _mm_shuffle_ps(p1, p2, _MM_SHUFFLE(2, 0, 2, 0));
312 auto ip = _mm_shuffle_ps(p1, p2, _MM_SHUFFLE(3, 1, 3, 1));
313
314 __m128 out;
315 f(rp, ip, out);
316
317 _mm_store_ps(output + i, out);
318 }
319 // deal with last partial packet
320 for (int i = n & (~3); i < n; ++i)
321 {
322 __m128 out;
323 f(_mm_set_ss(real(input[i])), _mm_set_ss(imag(input[i])), out);
324 output[i] = _mm_cvtss_f32(out);
325 }
326}

References SpareStrings::input, and SpareStrings::output.

◆ rotate_parallel_simd_aligned()

void simd_complex_conversions::rotate_parallel_simd_aligned ( const float *  oldPhase,
const float *  newPhase,
std::complex< float > *  output,
int  n 
)
inline

Definition at line 328 of file SimdComplexConversions_sse2.h.

331{
332 for (int i = 0; i <= n - 4; i += 4)
333 {
334 auto [sin, cos] = sincos_ps(
335 oldPhase ?
336 _mm_sub_ps(_mm_load_ps(newPhase + i), _mm_load_ps(oldPhase + i)) :
337 _mm_load_ps(newPhase + i));
338
339 // Safe according to C++ standard
340 auto p1 = _mm_load_ps(reinterpret_cast<float*>(output + i));
341 auto p2 = _mm_load_ps(reinterpret_cast<float*>(output + i + 2));
342
343 // p1 = {real(c1), imag(c1), real(c2), imag(c2)}
344 // p2 = {real(c3), imag(c3), real(c4), imag(c4)}
345
346 auto rp = _mm_shuffle_ps(p1, p2, _MM_SHUFFLE(2, 0, 2, 0));
347 auto ip = _mm_shuffle_ps(p1, p2, _MM_SHUFFLE(3, 1, 3, 1));
348
349
350 // We need to calculate (rp, ip) * (cos, sin) -> (rp*cos - ip*sin, rp*sin + ip*cos)
351
352 auto out_rp = _mm_sub_ps(_mm_mul_ps(rp, cos), _mm_mul_ps(ip, sin));
353 auto out_ip = _mm_add_ps(_mm_mul_ps(rp, sin), _mm_mul_ps(ip, cos));
354
355 p1 = _mm_unpacklo_ps(out_rp, out_ip);
356 p2 = _mm_unpackhi_ps(out_rp, out_ip);
357
358 _mm_store_ps(reinterpret_cast<float*>(output + i), p1);
359 _mm_store_ps(reinterpret_cast<float*>(output + i + 2), p2);
360 }
361 // deal with last partial packet
362 for (int i = n & (~3); i < n; ++i)
363 {
364 const auto theta = oldPhase ? newPhase[i] - oldPhase[i] : newPhase[i];
365 output[i] *= std::complex<float>(cosf(theta), sinf(theta));
366 }
367}
void sincos_ps(v4sf x, v4sf *s, v4sf *c)
Definition: SseMathFuncs.h:631

References SpareStrings::output, and sincos_ps().

Here is the call graph for this function:

◆ sincos_ps()

std::pair< __m128, __m128 > simd_complex_conversions::sincos_ps ( __m128  x)
inline

Definition at line 181 of file SimdComplexConversions_sse2.h.

182{
183 using namespace details;
184 __m128 xmm1, xmm2, xmm3 = _mm_setzero_ps(), sign_bit_sin, y;
185 __m128i emm0, emm2, emm4;
186
187 sign_bit_sin = x;
188 /* take the absolute value */
189 x = _mm_and_ps(x, _mm_set1_ps(inv_sign_mask));
190 /* extract the sign bit (upper one) */
191 sign_bit_sin = _mm_and_ps(sign_bit_sin, _mm_set1_ps(sign_mask));
192
193 /* scale by 4/Pi */
194 y = _mm_mul_ps(x, _mm_set1_ps(cephes_FOPI));
195
196 /* store the integer part of y in emm2 */
197 emm2 = _mm_cvttps_epi32(y);
198
199 /* j=(j+1) & (~1) (see the cephes sources) */
200 emm2 = _mm_add_epi32(emm2, _mm_set1_epi32(1));
201 emm2 = _mm_and_si128(emm2, _mm_set1_epi32(~1));
202 y = _mm_cvtepi32_ps(emm2);
203
204 emm4 = emm2;
205
206 /* get the swap sign flag for the sine */
207 emm0 = _mm_and_si128(emm2, _mm_set1_epi32(4));
208 emm0 = _mm_slli_epi32(emm0, 29);
209 __m128 swap_sign_bit_sin = _mm_castsi128_ps(emm0);
210
211 /* get the polynom selection mask for the sine*/
212 emm2 = _mm_and_si128(emm2, _mm_set1_epi32(2));
213 emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
214 __m128 poly_mask = _mm_castsi128_ps(emm2);
215
216 /* The magic pass: "Extended precision modular arithmetic"
217 x = ((x - y * DP1) - y * DP2) - y * DP3; */
218 xmm1 = _mm_set1_ps(minus_cephes_DP1);
219 xmm2 = _mm_set1_ps(minus_cephes_DP2);
220 xmm3 = _mm_set1_ps(minus_cephes_DP3);
221 xmm1 = _mm_mul_ps(y, xmm1);
222 xmm2 = _mm_mul_ps(y, xmm2);
223 xmm3 = _mm_mul_ps(y, xmm3);
224 x = _mm_add_ps(x, xmm1);
225 x = _mm_add_ps(x, xmm2);
226 x = _mm_add_ps(x, xmm3);
227
228 emm4 = _mm_sub_epi32(emm4, _mm_set1_epi32(2));
229 emm4 = _mm_andnot_si128(emm4, _mm_set1_epi32(4));
230 emm4 = _mm_slli_epi32(emm4, 29);
231 __m128 sign_bit_cos = _mm_castsi128_ps(emm4);
232
233 sign_bit_sin = _mm_xor_ps(sign_bit_sin, swap_sign_bit_sin);
234
235 /* Evaluate the first polynom (0 <= x <= Pi/4) */
236 __m128 z = _mm_mul_ps(x, x);
237 y = _mm_set1_ps(coscof_p0);
238
239 y = _mm_mul_ps(y, z);
240 y = _mm_add_ps(y, _mm_set1_ps(coscof_p1));
241 y = _mm_mul_ps(y, z);
242 y = _mm_add_ps(y, _mm_set1_ps(coscof_p2));
243 y = _mm_mul_ps(y, z);
244 y = _mm_mul_ps(y, z);
245 __m128 tmp = _mm_mul_ps(z, _mm_set1_ps(0.5f));
246 y = _mm_sub_ps(y, tmp);
247 y = _mm_add_ps(y, _mm_set1_ps(1));
248
249 /* Evaluate the second polynom (Pi/4 <= x <= 0) */
250
251 __m128 y2 = _mm_set1_ps(sincof_p0);
252 y2 = _mm_mul_ps(y2, z);
253 y2 = _mm_add_ps(y2, _mm_set1_ps(sincof_p1));
254 y2 = _mm_mul_ps(y2, z);
255 y2 = _mm_add_ps(y2, _mm_set1_ps(sincof_p2));
256 y2 = _mm_mul_ps(y2, z);
257 y2 = _mm_mul_ps(y2, x);
258 y2 = _mm_add_ps(y2, x);
259
260 /* select the correct result from the two polynoms */
261 xmm3 = poly_mask;
262 __m128 ysin2 = _mm_and_ps(xmm3, y2);
263 __m128 ysin1 = _mm_andnot_ps(xmm3, y);
264 y2 = _mm_sub_ps(y2, ysin2);
265 y = _mm_sub_ps(y, ysin1);
266
267 xmm1 = _mm_add_ps(ysin1, ysin2);
268 xmm2 = _mm_add_ps(y, y2);
269
270 /* update the sign */
271 return std::make_pair(
272 _mm_xor_ps(xmm1, sign_bit_sin), _mm_xor_ps(xmm2, sign_bit_cos));
273}

References simd_complex_conversions::details::cephes_FOPI, simd_complex_conversions::details::coscof_p0, simd_complex_conversions::details::coscof_p1, simd_complex_conversions::details::coscof_p2, simd_complex_conversions::details::inv_sign_mask, simd_complex_conversions::details::minus_cephes_DP1, simd_complex_conversions::details::minus_cephes_DP2, simd_complex_conversions::details::minus_cephes_DP3, simd_complex_conversions::details::sign_mask, simd_complex_conversions::details::sincof_p0, simd_complex_conversions::details::sincof_p1, and simd_complex_conversions::details::sincof_p2.

Referenced by rotate_parallel_simd_aligned(), and sincos_ss().

Here is the caller graph for this function:

◆ sincos_ss()

std::pair< float, float > simd_complex_conversions::sincos_ss ( float  angle)
inline

Definition at line 280 of file SimdComplexConversions_sse2.h.

281{
282 auto res = sincos_ps(_mm_set_ss(angle));
283 return std::make_pair(_mm_cvtss_f32(res.first), _mm_cvtss_f32(res.second));
284}

References sincos_ps().

Here is the call graph for this function:

◆ sqrt_ss()

float simd_complex_conversions::sqrt_ss ( float  x)
inline

Definition at line 291 of file SimdComplexConversions_sse2.h.

292{
293 __m128 sse_value = _mm_set_ss(x);
294 sse_value = _mm_sqrt_ss(sse_value);
295 return _mm_cvtss_f32(sse_value);
296}