Loading...
Searching...
No Matches
volk_common.h
Go to the documentation of this file.
/* -*- c++ -*- */
/*
* Copyright 2010, 2011, 2015-2017, 2019, 2020 Free Software Foundation, Inc.
* Copyright 2023 - 2025 Magnus Lundmark <magnuslundmark@gmail.com>
*
* This file is part of VOLK
*
* SPDX-License-Identifier: LGPL-3.0-or-later
*/
#ifndef INCLUDED_LIBVOLK_COMMON_H
#define INCLUDED_LIBVOLK_COMMON_H
// Cross-platform attribute macros
#if _MSC_VER
#define __VOLK_ATTR_ALIGNED(x) __declspec(align(x))
#define __VOLK_ATTR_UNUSED
#define __VOLK_ATTR_INLINE __forceinline
#define __VOLK_ATTR_DEPRECATED __declspec(deprecated)
#define __VOLK_ATTR_EXPORT __declspec(dllexport)
#define __VOLK_ATTR_IMPORT __declspec(dllimport)
#define __VOLK_PREFETCH(addr)
#define __VOLK_ASM __asm
#elif defined(__clang__)
// AppleClang also defines __GNUC__, so do this check first. These
// will probably be the same as for __GNUC__, but let's keep them
// separate just to be safe.
#define __VOLK_ATTR_ALIGNED(x) __attribute__((aligned(x)))
#define __VOLK_ATTR_UNUSED __attribute__((unused))
#define __VOLK_ATTR_INLINE __attribute__((always_inline))
#define __VOLK_ATTR_DEPRECATED __attribute__((deprecated))
#define __VOLK_ASM __asm__
#define __VOLK_ATTR_EXPORT __attribute__((visibility("default")))
#define __VOLK_ATTR_IMPORT __attribute__((visibility("default")))
#define __VOLK_PREFETCH(addr) __builtin_prefetch(addr)
#elif defined __GNUC__
#define __VOLK_ATTR_ALIGNED(x) __attribute__((aligned(x)))
#define __VOLK_ATTR_UNUSED __attribute__((unused))
#define __VOLK_ATTR_INLINE __attribute__((always_inline))
#define __VOLK_ATTR_DEPRECATED __attribute__((deprecated))
#define __VOLK_ASM __asm__
#if __GNUC__ >= 4
#define __VOLK_ATTR_EXPORT __attribute__((visibility("default")))
#define __VOLK_ATTR_IMPORT __attribute__((visibility("default")))
#else
#define __VOLK_ATTR_EXPORT
#define __VOLK_ATTR_IMPORT
#endif
#define __VOLK_PREFETCH(addr) __builtin_prefetch(addr)
#elif _MSC_VER
#define __VOLK_ATTR_ALIGNED(x) __declspec(align(x))
#define __VOLK_ATTR_UNUSED
#define __VOLK_ATTR_INLINE __forceinline
#define __VOLK_ATTR_DEPRECATED __declspec(deprecated)
#define __VOLK_ATTR_EXPORT __declspec(dllexport)
#define __VOLK_ATTR_IMPORT __declspec(dllimport)
#define __VOLK_PREFETCH(addr)
#define __VOLK_ASM __asm
#else
#define __VOLK_ATTR_ALIGNED(x)
#define __VOLK_ATTR_UNUSED
#define __VOLK_ATTR_INLINE
#define __VOLK_ATTR_DEPRECATED
#define __VOLK_ATTR_EXPORT
#define __VOLK_ATTR_IMPORT
#define __VOLK_PREFETCH(addr)
#define __VOLK_ASM __asm__
#endif
// Ignore annoying warnings in MSVC
#if defined(_MSC_VER)
#pragma warning(disable : 4244) //'conversion' conversion from 'type1' to 'type2',
// possible loss of data
#pragma warning(disable : 4305) //'identifier' : truncation from 'type1' to 'type2'
#endif
// C-linkage declaration macros
// FIXME: due to the usage of complex.h, require gcc for c-linkage
#if defined(__cplusplus) && (__GNUC__)
#define __VOLK_DECL_BEGIN extern "C" {
#define __VOLK_DECL_END }
#else
#define __VOLK_DECL_BEGIN
#define __VOLK_DECL_END
#endif
// Define VOLK_API for library symbols
// https://gcc.gnu.org/wiki/Visibility
#ifdef volk_EXPORTS
#define VOLK_API __VOLK_ATTR_EXPORT
#else
#define VOLK_API __VOLK_ATTR_IMPORT
#endif
// The bit128 union used by some
#include <stdint.h>
#ifdef LV_HAVE_SSE
#ifdef _WIN32
#include <intrin.h>
#else
#include <x86intrin.h>
#endif
#endif
union bit128 {
uint8_t i8[16];
uint16_t i16[8];
uint32_t i[4];
float f[4];
double d[2];
#ifdef LV_HAVE_SSE
__m128 float_vec;
#endif
#ifdef LV_HAVE_SSE2
__m128i int_vec;
__m128d double_vec;
#endif
};
union bit256 {
uint8_t i8[32];
uint16_t i16[16];
uint32_t i[8];
float f[8];
double d[4];
#ifdef LV_HAVE_AVX
__m256 float_vec;
__m256i int_vec;
__m256d double_vec;
#endif
};
#define bit128_p(x) ((union bit128*)(x))
#define bit256_p(x) ((union bit256*)(x))
// log2f
#include <math.h>
// +-Inf -> +-127.0f in order to match the behaviour of the SIMD kernels
// NaN -> NaN (preserved for consistency)
static inline float log2f_non_ieee(float f)
{
float const result = log2f(f);
// Return NaN for NaN inputs or negative values (preserves IEEE behavior for invalid
// inputs)
if (isnan(result))
return result;
// Map ±Inf to ±127.0f to match SIMD kernel behavior
return isinf(result) ? copysignf(127.0f, result) : result;
}
// Constant used to do log10 calculations as faster log2
// precalculated 10.0 / log2f_non_ieee(10.0) to allow for constexpr
#define volk_log2to10factor (0x1.815182p1) // 3.01029995663981209120
// arctan(x) polynomial expansion
static inline float volk_arctan_poly(const float x)
{
/*
* arctan(x) polynomial expansion on the interval [-1, 1]
* Maximum relative error < 6.6e-7
*/
const float a1 = +0x1.ffffeap-1f;
const float a3 = -0x1.55437p-2f;
const float a5 = +0x1.972be6p-3f;
const float a7 = -0x1.1436ap-3f;
const float a9 = +0x1.5785aap-4f;
const float a11 = -0x1.2f3004p-5f;
const float a13 = +0x1.01a37cp-7f;
const float x_times_x = x * x;
float arctan = a13;
arctan = fmaf(x_times_x, arctan, a11);
arctan = fmaf(x_times_x, arctan, a9);
arctan = fmaf(x_times_x, arctan, a7);
arctan = fmaf(x_times_x, arctan, a5);
arctan = fmaf(x_times_x, arctan, a3);
arctan = fmaf(x_times_x, arctan, a1);
arctan *= x;
return arctan;
}
// sin(x) polynomial expansion
static inline float volk_sin_poly(const float x)
{
/*
* Minimax polynomial for sin(x) on [-pi/4, pi/4]
* Coefficients via Remez algorithm (Sollya)
* Max |error| < 7.3e-9
* sin(x) = x + x^3 * (s1 + x^2 * (s2 + x^2 * s3))
*/
const float s1 = -0x1.555552p-3f;
const float s2 = +0x1.110be2p-7f;
const float s3 = -0x1.9ab22ap-13f;
const float x2 = x * x;
const float x3 = x2 * x;
float poly = fmaf(x2, s3, s2);
poly = fmaf(x2, poly, s1);
return fmaf(x3, poly, x);
}
// cos(x) polynomial expansion
static inline float volk_cos_poly(const float x)
{
/*
* Minimax polynomial for cos(x) on [-pi/4, pi/4]
* Coefficients via Remez algorithm (Sollya)
* Max |error| < 1.1e-7
* cos(x) = 1 + x^2 * (c1 + x^2 * (c2 + x^2 * c3))
*/
const float c1 = -0x1.fffff4p-2f;
const float c2 = +0x1.554a46p-5f;
const float c3 = -0x1.661be2p-10f;
const float x2 = x * x;
float poly = fmaf(x2, c3, c2);
poly = fmaf(x2, poly, c1);
return fmaf(x2, poly, 1.0f);
}
// sin(x) with Cody-Waite argument reduction
static inline float volk_sin(const float x)
{
/*
* Cody-Waite argument reduction: n = round(x * 2/pi), r = x - n * pi/2
* Then use sin/cos polynomials based on quadrant
*/
const float two_over_pi = 0x1.45f306p-1f;
const float pi_over_2_hi = 0x1.921fb6p+0f;
const float pi_over_2_lo = -0x1.777a5cp-25f;
float n_f = rintf(x * two_over_pi);
int n = (int)n_f;
float r = fmaf(-n_f, pi_over_2_hi, x);
r = fmaf(-n_f, pi_over_2_lo, r);
float sin_r = volk_sin_poly(r);
float cos_r = volk_cos_poly(r);
// Quadrant selection: n&1 swaps sin/cos, n&2 negates
float result = (n & 1) ? cos_r : sin_r;
return (n & 2) ? -result : result;
}
// cos(x) with Cody-Waite argument reduction
static inline float volk_cos(const float x)
{
/*
* Cody-Waite argument reduction: n = round(x * 2/pi), r = x - n * pi/2
* Then use sin/cos polynomials based on quadrant
*/
const float two_over_pi = 0x1.45f306p-1f;
const float pi_over_2_hi = 0x1.921fb6p+0f;
const float pi_over_2_lo = -0x1.777a5cp-25f;
float n_f = rintf(x * two_over_pi);
int n = (int)n_f;
float r = fmaf(-n_f, pi_over_2_hi, x);
r = fmaf(-n_f, pi_over_2_lo, r);
float sin_r = volk_sin_poly(r);
float cos_r = volk_cos_poly(r);
// Quadrant selection: n&1 swaps sin/cos, (n+1)&2 negates
float result = (n & 1) ? sin_r : cos_r;
return ((n + 1) & 2) ? -result : result;
}
// arctan(x)
static inline float volk_arctan(const float x)
{
/*
* arctan(x) + arctan(1 / x) == sign(x) * pi / 2
*/
const float pi_2 = 0x1.921fb6p0f;
// Propagate NaN
if (isnan(x)) {
return x;
}
// arctan(±∞) = ±π/2
if (isinf(x)) {
return copysignf(pi_2, x);
}
if (fabs(x) < 1.f) {
return volk_arctan_poly(x);
} else {
return copysignf(pi_2, x) - volk_arctan_poly(1.f / x);
}
}
// arctan2(y, x)
static inline float volk_atan2(const float y, const float x)
{
/*
* / arctan(y / x) if x > 0
* | arctan(y / x) + PI if x < 0 and y >= 0
* atan2(y, x) = | arctan(y / x) - PI if x < 0 and y < 0
* | sign(y) * PI / 2 if x = 0
* \ undefined if x = 0 and y = 0
* atan2f(0.f, 0.f) shall return 0.f
* atan2f(0.f, -0.f) shall return -0.f
*/
const float pi = 0x1.921fb6p1f;
const float pi_2 = 0x1.921fb6p0f;
// Propagate NaN from inputs
if (isnan(x) || isnan(y)) {
return x + y;
}
// Handle infinity cases per IEEE 754
if (isinf(y)) {
if (isinf(x)) {
// Both infinite: atan2(±∞, ±∞) = ±π/4 or ±3π/4
const float angle = (x > 0.f) ? (pi_2 / 2.f) : (3.f * pi_2 / 2.f);
return copysignf(angle, y);
} else {
// y infinite, x finite: atan2(±∞, x) = ±π/2
return copysignf(pi_2, y);
}
}
if (isinf(x)) {
// x infinite, y finite: atan2(y, +∞) = ±0, atan2(y, -∞) = ±π
return (x > 0.f) ? copysignf(0.f, y) : copysignf(pi, y);
}
if (fabs(x) == 0.f) {
return (fabs(y) == 0.f) ? copysignf(0.f, y) : copysignf(pi_2, y);
}
const int swap = fabs(x) < fabs(y);
const float numerator = swap ? x : y;
const float denominator = swap ? y : x;
float input = numerator / denominator;
if (isnan(input)) {
input = numerator;
}
float result = volk_arctan_poly(input);
result = swap ? (input >= 0.f ? pi_2 : -pi_2) - result : result;
if (x < 0.f) {
result += copysignf(pi, y);
}
return result;
}
// arcsin(x) polynomial expansion
// P(u) such that asin(x) = x * P(x^2) on |x| <= 0.5
// Maximum relative error ~1.5e-6
static inline float volk_arcsin_poly(const float x)
{
const float c0 = 0x1.ffffcep-1f;
const float c1 = 0x1.55b648p-3f;
const float c2 = 0x1.24d192p-4f;
const float c3 = 0x1.0a788p-4f;
const float u = x * x;
float p = c3;
p = fmaf(u, p, c2);
p = fmaf(u, p, c1);
p = fmaf(u, p, c0);
return x * p;
}
// arcsin(x) using two-range algorithm
static inline float volk_arcsin(const float x)
{
const float pi_2 = 0x1.921fb6p0f;
const float ax = fabsf(x);
if (ax <= 0.5f) {
// Small argument: direct polynomial
return volk_arcsin_poly(x);
} else {
// Large argument: use identity asin(x) = pi/2 - 2*asin(sqrt((1-|x|)/2))
const float t = (1.0f - ax) * 0.5f;
const float s = sqrtf(t);
const float inner = volk_arcsin_poly(s);
const float result = pi_2 - 2.0f * inner;
return copysignf(result, x);
}
}
// arccos(x) = pi/2 - arcsin(x)
static inline float volk_arccos(const float x)
{
const float pi_2 = 0x1.921fb6p0f;
return pi_2 - volk_arcsin(x);
}
#endif /*INCLUDED_LIBVOLK_COMMON_H*/