Loading...
Searching...
No Matches
volk_32f_asin_32f.h
Go to the documentation of this file.
/* -*- c++ -*- */
/*
* Copyright 2014 Free Software Foundation, Inc.
* Copyright 2025 Magnus Lundmark <magnuslundmark@gmail.com>
*
* This file is part of VOLK
*
* SPDX-License-Identifier: LGPL-3.0-or-later
*/
#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include <volk/volk_common.h>
#ifndef INCLUDED_volk_32f_asin_32f_a_H
#define INCLUDED_volk_32f_asin_32f_a_H
#ifdef LV_HAVE_GENERIC
static inline void
volk_32f_asin_32f_generic(float* bVector, const float* aVector, unsigned int num_points)
{
for (unsigned int i = 0; i < num_points; i++) {
bVector[i] = volk_arcsin(aVector[i]);
}
}
#endif /* LV_HAVE_GENERIC */
#ifdef LV_HAVE_SSE4_1
#include <smmintrin.h>
#include <volk/volk_sse_intrinsics.h>
static inline void
volk_32f_asin_32f_a_sse4_1(float* bVector, const float* aVector, unsigned int num_points)
{
const __m128 pi_2 = _mm_set1_ps(0x1.921fb6p0f);
const __m128 half = _mm_set1_ps(0.5f);
const __m128 one = _mm_set1_ps(1.0f);
const __m128 two = _mm_set1_ps(2.0f);
const __m128 sign_mask = _mm_set1_ps(-0.0f);
unsigned int number = 0;
const unsigned int quarterPoints = num_points / 4;
for (; number < quarterPoints; number++) {
__m128 aVal = _mm_load_ps(aVector);
// Get absolute value and sign
__m128 sign = _mm_and_ps(aVal, sign_mask);
__m128 ax = _mm_andnot_ps(sign_mask, aVal);
// Two-range computation
// Small: result = arcsin_poly(x)
// Large: result = pi/2 - 2*arcsin_poly(sqrt((1-|x|)/2))
__m128 t = _mm_mul_ps(_mm_sub_ps(one, ax), half);
__m128 s = _mm_sqrt_ps(t);
// Compute polynomial for both ranges
__m128 poly_small = _mm_arcsin_poly_sse(ax);
__m128 poly_large = _mm_arcsin_poly_sse(s);
// Large range: pi/2 - 2*poly_large
__m128 result_large = _mm_sub_ps(pi_2, _mm_mul_ps(two, poly_large));
// Blend based on |x| > 0.5
__m128 mask = _mm_cmpgt_ps(ax, half);
__m128 result = _mm_blendv_ps(poly_small, result_large, mask);
// Apply sign
result = _mm_or_ps(result, sign);
_mm_store_ps(bVector, result);
aVector += 4;
bVector += 4;
}
number = quarterPoints * 4;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_SSE4_1 */
#ifdef LV_HAVE_AVX
#include <immintrin.h>
#include <volk/volk_avx_intrinsics.h>
static inline void
volk_32f_asin_32f_a_avx(float* bVector, const float* aVector, unsigned int num_points)
{
const __m256 pi_2 = _mm256_set1_ps(0x1.921fb6p0f);
const __m256 half = _mm256_set1_ps(0.5f);
const __m256 one = _mm256_set1_ps(1.0f);
const __m256 two = _mm256_set1_ps(2.0f);
const __m256 sign_mask = _mm256_set1_ps(-0.0f);
unsigned int number = 0;
const unsigned int eighthPoints = num_points / 8;
for (; number < eighthPoints; number++) {
__m256 aVal = _mm256_load_ps(aVector);
// Get absolute value and sign
__m256 sign = _mm256_and_ps(aVal, sign_mask);
__m256 ax = _mm256_andnot_ps(sign_mask, aVal);
// Two-range computation
__m256 t = _mm256_mul_ps(_mm256_sub_ps(one, ax), half);
__m256 s = _mm256_sqrt_ps(t);
// Compute polynomial for both ranges
__m256 poly_small = _mm256_arcsin_poly_avx(ax);
__m256 poly_large = _mm256_arcsin_poly_avx(s);
// Large range: pi/2 - 2*poly_large
__m256 result_large = _mm256_sub_ps(pi_2, _mm256_mul_ps(two, poly_large));
// Blend based on |x| > 0.5
__m256 mask = _mm256_cmp_ps(ax, half, _CMP_GT_OS);
__m256 result = _mm256_blendv_ps(poly_small, result_large, mask);
// Apply sign
result = _mm256_or_ps(result, sign);
_mm256_store_ps(bVector, result);
aVector += 8;
bVector += 8;
}
number = eighthPoints * 8;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_AVX */
#ifdef LV_HAVE_AVX2
#include <immintrin.h>
#include <volk/volk_avx2_fma_intrinsics.h>
static inline void volk_32f_asin_32f_a_avx2_fma(float* bVector,
const float* aVector,
unsigned int num_points)
{
const __m256 pi_2 = _mm256_set1_ps(0x1.921fb6p0f);
const __m256 half = _mm256_set1_ps(0.5f);
const __m256 one = _mm256_set1_ps(1.0f);
const __m256 two = _mm256_set1_ps(2.0f);
const __m256 sign_mask = _mm256_set1_ps(-0.0f);
unsigned int number = 0;
const unsigned int eighthPoints = num_points / 8;
for (; number < eighthPoints; number++) {
__m256 aVal = _mm256_load_ps(aVector);
// Get absolute value and sign
__m256 sign = _mm256_and_ps(aVal, sign_mask);
__m256 ax = _mm256_andnot_ps(sign_mask, aVal);
// Two-range computation
__m256 t = _mm256_mul_ps(_mm256_sub_ps(one, ax), half);
__m256 s = _mm256_sqrt_ps(t);
// Compute polynomial for both ranges
__m256 poly_small = _mm256_arcsin_poly_avx2_fma(ax);
__m256 poly_large = _mm256_arcsin_poly_avx2_fma(s);
// Large range: pi/2 - 2*poly_large
__m256 result_large = _mm256_fnmadd_ps(two, poly_large, pi_2);
// Blend based on |x| > 0.5
__m256 mask = _mm256_cmp_ps(ax, half, _CMP_GT_OS);
__m256 result = _mm256_blendv_ps(poly_small, result_large, mask);
// Apply sign
result = _mm256_or_ps(result, sign);
_mm256_store_ps(bVector, result);
aVector += 8;
bVector += 8;
}
number = eighthPoints * 8;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_AVX2 */
#ifdef LV_HAVE_AVX512F
#include <immintrin.h>
#include <volk/volk_avx512_intrinsics.h>
static inline void
volk_32f_asin_32f_a_avx512(float* bVector, const float* aVector, unsigned int num_points)
{
const __m512 pi_2 = _mm512_set1_ps(0x1.921fb6p0f);
const __m512 half = _mm512_set1_ps(0.5f);
const __m512 one = _mm512_set1_ps(1.0f);
const __m512 two = _mm512_set1_ps(2.0f);
const __m512i sign_mask = _mm512_set1_epi32(0x80000000);
unsigned int number = 0;
const unsigned int sixteenthPoints = num_points / 16;
for (; number < sixteenthPoints; number++) {
__m512 aVal = _mm512_load_ps(aVector);
// Get absolute value and sign using integer ops (AVX512F compatible)
__m512i aVal_i = _mm512_castps_si512(aVal);
__m512i sign = _mm512_and_epi32(aVal_i, sign_mask);
__m512 ax = _mm512_castsi512_ps(_mm512_andnot_epi32(sign_mask, aVal_i));
// Two-range computation
__m512 t = _mm512_mul_ps(_mm512_sub_ps(one, ax), half);
__m512 s = _mm512_sqrt_ps(t);
// Compute polynomial for both ranges
__m512 poly_small = _mm512_arcsin_poly_avx512(ax);
__m512 poly_large = _mm512_arcsin_poly_avx512(s);
// Large range: pi/2 - 2*poly_large
__m512 result_large = _mm512_fnmadd_ps(two, poly_large, pi_2);
// Blend based on |x| > 0.5
__mmask16 mask = _mm512_cmp_ps_mask(ax, half, _CMP_GT_OS);
__m512 result = _mm512_mask_blend_ps(mask, poly_small, result_large);
// Apply sign
result = _mm512_castsi512_ps(_mm512_or_epi32(_mm512_castps_si512(result), sign));
_mm512_store_ps(bVector, result);
aVector += 16;
bVector += 16;
}
number = sixteenthPoints * 16;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_AVX512F */
#endif /* INCLUDED_volk_32f_asin_32f_a_H */
#ifndef INCLUDED_volk_32f_asin_32f_u_H
#define INCLUDED_volk_32f_asin_32f_u_H
#ifdef LV_HAVE_SSE4_1
#include <smmintrin.h>
#include <volk/volk_sse_intrinsics.h>
static inline void
volk_32f_asin_32f_u_sse4_1(float* bVector, const float* aVector, unsigned int num_points)
{
const __m128 pi_2 = _mm_set1_ps(0x1.921fb6p0f);
const __m128 half = _mm_set1_ps(0.5f);
const __m128 one = _mm_set1_ps(1.0f);
const __m128 two = _mm_set1_ps(2.0f);
const __m128 sign_mask = _mm_set1_ps(-0.0f);
unsigned int number = 0;
const unsigned int quarterPoints = num_points / 4;
for (; number < quarterPoints; number++) {
__m128 aVal = _mm_loadu_ps(aVector);
__m128 sign = _mm_and_ps(aVal, sign_mask);
__m128 ax = _mm_andnot_ps(sign_mask, aVal);
__m128 t = _mm_mul_ps(_mm_sub_ps(one, ax), half);
__m128 s = _mm_sqrt_ps(t);
__m128 poly_small = _mm_arcsin_poly_sse(ax);
__m128 poly_large = _mm_arcsin_poly_sse(s);
__m128 result_large = _mm_sub_ps(pi_2, _mm_mul_ps(two, poly_large));
__m128 mask = _mm_cmpgt_ps(ax, half);
__m128 result = _mm_blendv_ps(poly_small, result_large, mask);
result = _mm_or_ps(result, sign);
_mm_storeu_ps(bVector, result);
aVector += 4;
bVector += 4;
}
number = quarterPoints * 4;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_SSE4_1 */
#ifdef LV_HAVE_AVX
#include <immintrin.h>
#include <volk/volk_avx_intrinsics.h>
static inline void
volk_32f_asin_32f_u_avx(float* bVector, const float* aVector, unsigned int num_points)
{
const __m256 pi_2 = _mm256_set1_ps(0x1.921fb6p0f);
const __m256 half = _mm256_set1_ps(0.5f);
const __m256 one = _mm256_set1_ps(1.0f);
const __m256 two = _mm256_set1_ps(2.0f);
const __m256 sign_mask = _mm256_set1_ps(-0.0f);
unsigned int number = 0;
const unsigned int eighthPoints = num_points / 8;
for (; number < eighthPoints; number++) {
__m256 aVal = _mm256_loadu_ps(aVector);
__m256 sign = _mm256_and_ps(aVal, sign_mask);
__m256 ax = _mm256_andnot_ps(sign_mask, aVal);
__m256 t = _mm256_mul_ps(_mm256_sub_ps(one, ax), half);
__m256 s = _mm256_sqrt_ps(t);
__m256 poly_small = _mm256_arcsin_poly_avx(ax);
__m256 poly_large = _mm256_arcsin_poly_avx(s);
__m256 result_large = _mm256_sub_ps(pi_2, _mm256_mul_ps(two, poly_large));
__m256 mask = _mm256_cmp_ps(ax, half, _CMP_GT_OS);
__m256 result = _mm256_blendv_ps(poly_small, result_large, mask);
result = _mm256_or_ps(result, sign);
_mm256_storeu_ps(bVector, result);
aVector += 8;
bVector += 8;
}
number = eighthPoints * 8;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_AVX */
#ifdef LV_HAVE_AVX2
#include <immintrin.h>
#include <volk/volk_avx2_fma_intrinsics.h>
static inline void volk_32f_asin_32f_u_avx2_fma(float* bVector,
const float* aVector,
unsigned int num_points)
{
const __m256 pi_2 = _mm256_set1_ps(0x1.921fb6p0f);
const __m256 half = _mm256_set1_ps(0.5f);
const __m256 one = _mm256_set1_ps(1.0f);
const __m256 two = _mm256_set1_ps(2.0f);
const __m256 sign_mask = _mm256_set1_ps(-0.0f);
unsigned int number = 0;
const unsigned int eighthPoints = num_points / 8;
for (; number < eighthPoints; number++) {
__m256 aVal = _mm256_loadu_ps(aVector);
__m256 sign = _mm256_and_ps(aVal, sign_mask);
__m256 ax = _mm256_andnot_ps(sign_mask, aVal);
__m256 t = _mm256_mul_ps(_mm256_sub_ps(one, ax), half);
__m256 s = _mm256_sqrt_ps(t);
__m256 poly_small = _mm256_arcsin_poly_avx2_fma(ax);
__m256 poly_large = _mm256_arcsin_poly_avx2_fma(s);
__m256 result_large = _mm256_fnmadd_ps(two, poly_large, pi_2);
__m256 mask = _mm256_cmp_ps(ax, half, _CMP_GT_OS);
__m256 result = _mm256_blendv_ps(poly_small, result_large, mask);
result = _mm256_or_ps(result, sign);
_mm256_storeu_ps(bVector, result);
aVector += 8;
bVector += 8;
}
number = eighthPoints * 8;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_AVX2 */
#ifdef LV_HAVE_AVX512F
#include <immintrin.h>
#include <volk/volk_avx512_intrinsics.h>
static inline void
volk_32f_asin_32f_u_avx512(float* bVector, const float* aVector, unsigned int num_points)
{
const __m512 pi_2 = _mm512_set1_ps(0x1.921fb6p0f);
const __m512 half = _mm512_set1_ps(0.5f);
const __m512 one = _mm512_set1_ps(1.0f);
const __m512 two = _mm512_set1_ps(2.0f);
const __m512i sign_mask = _mm512_set1_epi32(0x80000000);
unsigned int number = 0;
const unsigned int sixteenthPoints = num_points / 16;
for (; number < sixteenthPoints; number++) {
__m512 aVal = _mm512_loadu_ps(aVector);
__m512i aVal_i = _mm512_castps_si512(aVal);
__m512i sign = _mm512_and_epi32(aVal_i, sign_mask);
__m512 ax = _mm512_castsi512_ps(_mm512_andnot_epi32(sign_mask, aVal_i));
__m512 t = _mm512_mul_ps(_mm512_sub_ps(one, ax), half);
__m512 s = _mm512_sqrt_ps(t);
__m512 poly_small = _mm512_arcsin_poly_avx512(ax);
__m512 poly_large = _mm512_arcsin_poly_avx512(s);
__m512 result_large = _mm512_fnmadd_ps(two, poly_large, pi_2);
__mmask16 mask = _mm512_cmp_ps_mask(ax, half, _CMP_GT_OS);
__m512 result = _mm512_mask_blend_ps(mask, poly_small, result_large);
result = _mm512_castsi512_ps(_mm512_or_epi32(_mm512_castps_si512(result), sign));
_mm512_storeu_ps(bVector, result);
aVector += 16;
bVector += 16;
}
number = sixteenthPoints * 16;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_AVX512F */
#ifdef LV_HAVE_NEON
#include <arm_neon.h>
#include <volk/volk_neon_intrinsics.h>
static inline void
volk_32f_asin_32f_neon(float* bVector, const float* aVector, unsigned int num_points)
{
const float32x4_t pi_2 = vdupq_n_f32(0x1.921fb6p0f);
const float32x4_t half = vdupq_n_f32(0.5f);
const float32x4_t one = vdupq_n_f32(1.0f);
const float32x4_t two = vdupq_n_f32(2.0f);
unsigned int number = 0;
const unsigned int quarterPoints = num_points / 4;
for (; number < quarterPoints; number++) {
float32x4_t aVal = vld1q_f32(aVector);
// Get absolute value and sign
float32x4_t ax = vabsq_f32(aVal);
uint32x4_t sign_bits =
vandq_u32(vreinterpretq_u32_f32(aVal), vdupq_n_u32(0x80000000));
// Two-range computation
float32x4_t t = vmulq_f32(vsubq_f32(one, ax), half);
float32x4_t s = _vsqrtq_f32(t);
// Compute polynomial for both ranges
float32x4_t poly_small = _varcsinq_f32(ax);
float32x4_t poly_large = _varcsinq_f32(s);
// Large range: pi/2 - 2*poly_large
float32x4_t result_large = vmlsq_f32(pi_2, two, poly_large);
// Blend based on |x| > 0.5
uint32x4_t mask = vcgtq_f32(ax, half);
float32x4_t result = vbslq_f32(mask, result_large, poly_small);
// Apply sign
result =
vreinterpretq_f32_u32(vorrq_u32(vreinterpretq_u32_f32(result), sign_bits));
vst1q_f32(bVector, result);
aVector += 4;
bVector += 4;
}
number = quarterPoints * 4;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_NEON */
#ifdef LV_HAVE_NEONV8
#include <arm_neon.h>
#include <volk/volk_neon_intrinsics.h>
static inline void
volk_32f_asin_32f_neonv8(float* bVector, const float* aVector, unsigned int num_points)
{
const float32x4_t pi_2 = vdupq_n_f32(0x1.921fb6p0f);
const float32x4_t half = vdupq_n_f32(0.5f);
const float32x4_t one = vdupq_n_f32(1.0f);
const float32x4_t two = vdupq_n_f32(2.0f);
unsigned int number = 0;
const unsigned int quarterPoints = num_points / 4;
for (; number < quarterPoints; number++) {
float32x4_t aVal = vld1q_f32(aVector);
float32x4_t ax = vabsq_f32(aVal);
uint32x4_t sign_bits =
vandq_u32(vreinterpretq_u32_f32(aVal), vdupq_n_u32(0x80000000));
float32x4_t t = vmulq_f32(vsubq_f32(one, ax), half);
float32x4_t s = vsqrtq_f32(t);
float32x4_t poly_small = _varcsinq_f32_neonv8(ax);
float32x4_t poly_large = _varcsinq_f32_neonv8(s);
float32x4_t result_large = vfmsq_f32(pi_2, two, poly_large);
uint32x4_t mask = vcgtq_f32(ax, half);
float32x4_t result = vbslq_f32(mask, result_large, poly_small);
result =
vreinterpretq_f32_u32(vorrq_u32(vreinterpretq_u32_f32(result), sign_bits));
vst1q_f32(bVector, result);
aVector += 4;
bVector += 4;
}
number = quarterPoints * 4;
for (; number < num_points; number++) {
*bVector++ = volk_arcsin(*aVector++);
}
}
#endif /* LV_HAVE_NEONV8 */
#endif /* INCLUDED_volk_32f_asin_32f_u_H */