Switched to abseil for span (they're more verbose but disappear better with compilers)

This commit is contained in:
paul 2019-08-21 18:18:12 +02:00
parent 31cb328969
commit 2a9991dfaa
14 changed files with 1368 additions and 107 deletions

View file

@ -140,7 +140,7 @@ add_executable(sfizz_tests ${TEST_SOURCES} ${COMMON_SOURCES})
if(UNIX)
target_link_libraries(sfizz_tests stdc++fs)
endif(UNIX)
target_link_libraries(sfizz_tests Catch2::Catch2 absl::strings absl::flat_hash_map sndfile readerwriterqueue cnpy-static gsl::gsl-lite)
target_link_libraries(sfizz_tests Catch2::Catch2 absl::strings absl::flat_hash_map sndfile readerwriterqueue cnpy-static absl::span)
target_include_directories(sfizz_tests SYSTEM PRIVATE sources)
file(COPY "tests" DESTINATION ${CMAKE_BINARY_DIR})
@ -150,7 +150,7 @@ target_link_libraries(bench_span benchmark gsl::gsl-lite absl::span)
###############################
add_executable(bench_opf_high_vs_low benchmarks/OPF_high_vs_low.cpp)
target_link_libraries(bench_opf_high_vs_low benchmark gsl::gsl-lite)
target_link_libraries(bench_opf_high_vs_low benchmark absl::span)
###############################
add_executable(bench_looping_index benchmarks/Looping_index.cpp)
@ -166,7 +166,19 @@ add_executable(bm_cum_sum benchmarks/Cum_Sum.cpp)
target_link_libraries(bm_cum_sum benchmark)
add_executable(bm_write benchmarks/BM_writeInterleaved.cpp sources/SIMDSSE.cpp)
target_link_libraries(bm_write benchmark gsl::gsl-lite)
target_link_libraries(bm_write benchmark absl::span)
add_executable(bm_fill benchmarks/BM_fill.cpp sources/SIMDSSE.cpp)
target_link_libraries(bm_fill benchmark gsl::gsl-lite)
target_link_libraries(bm_fill benchmark absl::span)
add_executable(bm_math_functions benchmarks/Math_functions.cpp)
target_link_libraries(bm_math_functions benchmark)
if (UNIX)
target_compile_options(bm_math_functions PRIVATE -funsafe-math-optimizations -msse4.2)
endif()
add_executable(bm_mathfuns benchmarks/BM_mathfuns.cpp sources/SIMDSSE.cpp)
target_link_libraries(bm_mathfuns benchmark absl::span absl::span)
if (UNIX)
# target_compile_options(bm_math_loops PRIVATE -fopenmp)
endif()

View file

@ -0,0 +1,87 @@
#include <benchmark/benchmark.h>
#include <random>
#include <numeric>
#include <vector>
#include <cmath>
#include <iostream>
#include "../sources/SIMDHelpers.h"
static void Dummy(benchmark::State& state) {
std::random_device rd { };
std::mt19937 gen { rd() };
std::vector<float> source(state.range(0));
std::vector<float> result(state.range(0));
std::normal_distribution<float> dist { };
std::generate(source.begin(), source.end(), [&]() { return dist(gen); });
for (auto _ : state)
{
for (int i = 0; i < state.range(0); ++i)
result[i] = source[i];
benchmark::DoNotOptimize(result);
}
}
static void StdExp(benchmark::State& state) {
std::random_device rd { };
std::mt19937 gen { rd() };
std::vector<float> source(state.range(0));
std::vector<float> result(state.range(0));
std::normal_distribution<float> dist { };
std::generate(source.begin(), source.end(), [&]() { return dist(gen); });
for (auto _ : state)
{
for (int i = 0; i < state.range(0); ++i)
result[i] = std::exp(source[i]);
benchmark::DoNotOptimize(result);
}
}
// static void StdExpOMP(benchmark::State& state) {
// std::random_device rd { };
// std::mt19937 gen { rd() };
// std::vector<float> source(state.range(0));
// std::vector<float> result(state.range(0));
// std::normal_distribution<float> dist { };
// std::generate(source.begin(), source.end(), [&]() { return dist(gen); });
// for (auto _ : state)
// {
// #pragma omp simd
// for (int i = 0; i < state.range(0); ++i)
// result[i] = std::exp(source[i]);
// benchmark::DoNotOptimize(result);
// }
// }
static void Scalar(benchmark::State& state) {
std::random_device rd { };
std::mt19937 gen { rd() };
std::vector<float> source(state.range(0));
std::vector<float> result(state.range(0));
std::normal_distribution<float> dist { };
std::generate(source.begin(), source.end(), [&]() { return dist(gen); });
for (auto _ : state)
{
exp<float, false>(source, absl::MakeSpan(result));
benchmark::DoNotOptimize(result);
}
}
static void SIMD(benchmark::State& state) {
std::random_device rd { };
std::mt19937 gen { rd() };
std::vector<float> source(state.range(0));
std::vector<float> result(state.range(0));
std::normal_distribution<float> dist { };
std::generate(source.begin(), source.end(), [&]() { return dist(gen); });
for (auto _ : state)
{
exp<float, true>(source, absl::MakeSpan(result));
benchmark::DoNotOptimize(result);
}
}
BENCHMARK(Dummy)->RangeMultiplier(2)->Range(1 << 6, 1 << 10);
BENCHMARK(StdExp)->RangeMultiplier(2)->Range(1 << 6, 1 << 10);
BENCHMARK(Scalar)->RangeMultiplier(2)->Range(1 << 6, 1 << 10);
BENCHMARK(SIMD)->RangeMultiplier(2)->Range(1 << 6, 1 << 10);
BENCHMARK_MAIN();

View file

@ -0,0 +1,111 @@
#include <benchmark/benchmark.h>
#include <random>
#include <cmath>
#include <iostream>
/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1984, 1987, 1989 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Single precision exponential function.
* test interval: [-0.5, +0.5]
* trials: 80000
* peak relative error: 7.6e-8
* rms relative error: 2.8e-8
*/
static float MAXNUMF = 3.4028234663852885981170418348451692544e38;
static float MAXLOGF = 88.72283905206835;
static float MINLOGF = -103.278929903431851103; /* log(2^-149) */
static float LOG2EF = 1.44269504088896341;
static float C1 = 0.693359375;
static float C2 = -2.12194440e-4;
float cephes_expf(float xx) {
float x, z;
int n;
x = xx;
if( x > MAXLOGF)
{
//mtherr( "expf", OVERFLOW );
return( MAXNUMF );
}
if( x < MINLOGF )
{
//mtherr( "expf", UNDERFLOW );
return(0.0);
}
/* Express e**x = e**g 2**n
* = e**g e**( n loge(2) )
* = e**( g + n loge(2) )
*/
z = floorf( LOG2EF * x + 0.5 ); /* floor() truncates toward -infinity. */
x -= z * C1;
x -= z * C2;
n = z;
z = x * x;
/* Theoretical peak relative error in [-0.5, +0.5] is 4.2e-9. */
z =
((((( 1.9875691500E-4f * x
+ 1.3981999507E-3f) * x
+ 8.3334519073E-3f) * x
+ 4.1665795894E-2f) * x
+ 1.6666665459E-1f) * x
+ 5.0000001201E-1f) * z
+ x
+ 1.0;
/* multiply by power of 2 */
x = ldexpf( z, n );
return( x );
}
static void Dummy(benchmark::State& state) {
std::random_device rd { };
std::mt19937 gen { rd() };
std::normal_distribution<float> dist { };
for (auto _ : state)
{
auto value = dist(gen);
benchmark::DoNotOptimize(value);
}
}
static void StdExp(benchmark::State& state) {
std::random_device rd { };
std::mt19937 gen { rd() };
std::normal_distribution<float> dist { };
for (auto _ : state)
{
auto value = std::exp(dist(gen));
benchmark::DoNotOptimize(value);
}
}
static void CephesExp(benchmark::State& state) {
std::random_device rd { };
std::mt19937 gen { rd() };
std::normal_distribution<float> dist { };
for (auto _ : state)
{
auto value = cephes_expf(dist(gen));
benchmark::DoNotOptimize(value);
}
}
BENCHMARK(Dummy);
BENCHMARK(StdExp);
BENCHMARK(CephesExp);
BENCHMARK_MAIN();

View file

@ -1,5 +1,5 @@
#include "FilePool.h"
#include "gsl/gsl-lite.hpp"
#include <chrono>
using namespace std::chrono_literals;
@ -24,7 +24,7 @@ std::optional<sfz::FilePool::FileInformation> sfz::FilePool::getFileInformation(
auto preloadedSize = std::min(returnedValue.end, static_cast<uint32_t>(config::preloadSize));
returnedValue.preloadedData = std::make_shared<StereoBuffer<float>>(preloadedSize);
sndFile.readf(tempReadBuffer.data(), preloadedSize);
returnedValue.preloadedData->readInterleaved(gsl::make_span(tempReadBuffer).first(preloadedSize));
returnedValue.preloadedData->readInterleaved(absl::MakeSpan(tempReadBuffer).first(preloadedSize));
preloadedData[filename] = returnedValue.preloadedData;
// char buffer [2048] ;
// sndFile.command(SFC_GET_LOG_INFO, buffer, sizeof(buffer)) ;

View file

@ -1,6 +1,6 @@
#include "Globals.h"
#include <cmath>
#include "gsl/gsl-lite.hpp"
#include <absl/types/span.h>
template<class Type=float>
class OnePoleFilter
@ -27,7 +27,7 @@ public:
Type getGain() const { return gain; }
int processLowpass(gsl::span<const Type> input, gsl::span<Type> lowpass)
int processLowpass(absl::Span<const Type> input, absl::Span<Type> lowpass)
{
for (auto [in, out] = std::pair(input.begin(), lowpass.begin());
in < input.end() && out < lowpass.end(); in++, out++)
@ -37,7 +37,7 @@ public:
return std::min(input.size(), lowpass.size());
}
int processHighpass(gsl::span<const Type> input, gsl::span<Type> highpass)
int processHighpass(absl::Span<const Type> input, absl::Span<Type> highpass)
{
for (auto [in, out] = std::pair(input.begin(), highpass.begin());
in < input.end() && out < highpass.end(); in++, out++)
@ -47,7 +47,7 @@ public:
return std::min(input.size(), highpass.size());
}
int processLowpassVariableGain(gsl::span<const Type> input, gsl::span<Type> lowpass, gsl::span<const Type> gain)
int processLowpassVariableGain(absl::Span<const Type> input, absl::Span<Type> lowpass, absl::Span<const Type> gain)
{
for (auto [in, out, g] = std::tuple(input.begin(), lowpass.begin(), gain.begin());
in < input.end() && out < lowpass.end() && g < gain.end(); in++, out++, g++)
@ -59,7 +59,7 @@ public:
return std::min({ input.size(), lowpass.size(), gain.size() });
}
int processHighpassVariableGain(gsl::span<const Type> input, gsl::span<Type> highpass, gsl::span<const Type> gain)
int processHighpassVariableGain(absl::Span<const Type> input, absl::Span<Type> highpass, absl::Span<const Type> gain)
{
for (auto [in, out, g] = std::tuple(input.begin(), highpass.begin(), gain.begin());
in < input.end() && out < highpass.end() && g < gain.end(); in++, out++, g++)

View file

@ -2,34 +2,40 @@
#include "Helpers.h"
template<>
void readInterleaved<float, true>(gsl::span<const float> input, gsl::span<float> outputLeft, gsl::span<float> outputRight) noexcept
void readInterleaved<float, true>(absl::Span<const float> input, absl::Span<float> outputLeft, absl::Span<float> outputRight) noexcept
{
readInterleaved<float, false>(input, outputLeft, outputRight);
}
template<>
void writeInterleaved<float, true>(gsl::span<const float> inputLeft, gsl::span<const float> inputRight, gsl::span<float> output) noexcept
void writeInterleaved<float, true>(absl::Span<const float> inputLeft, absl::Span<const float> inputRight, absl::Span<float> output) noexcept
{
writeInterleaved<float, false>(inputLeft, inputRight, output);
}
// template<class Type, bool SIMD=false>
// void loopingSFZIndex(gsl::span<const Type> inputLeft, gsl::span<const Type> inputRight, gsl::span<Type> output);
// void loopingSFZIndex(absl::Span<const Type> inputLeft, absl::Span<const Type> inputRight, absl::Span<Type> output);
// template<class Type, bool SIMD=false>
// void linearRamp(gsl::span<Type> output, Type start, Type end);
// void linearRamp(absl::Span<Type> output, Type start, Type end);
// template<class Type, bool SIMD=false>
// void exponentialRamp(gsl::span<Type> output, Type start, Type end);
// void exponentialRamp(absl::Span<Type> output, Type start, Type end);
// template<class Type, bool SIMD=false>
// void applyGain(Type gain, gsl::span<Type> output);
// void applyGain(Type gain, absl::Span<Type> output);
// template<class Type, bool SIMD=false>
// void applyGain(gsl::span<const Type> output, gsl::span<Type> output);
// void applyGain(absl::Span<const Type> output, absl::Span<Type> output);
template<>
void fill<float, true>(gsl::span<float> output, float value) noexcept
void fill<float, true>(absl::Span<float> output, float value) noexcept
{
fill<float, false>(output, value);
}
template<>
void exp<float, true>(absl::Span<const float> input, absl::Span<float> output) noexcept
{
exp<float, false>(input, output);
}

View file

@ -1,9 +1,11 @@
#include "gsl/gsl-lite.hpp"
#include "absl/types/span.h"
#include "Globals.h"
#include "Helpers.h"
#include <cmath>
template<class T, bool SIMD=SIMDConfig::readInterleaved>
void readInterleaved(gsl::span<const T> input, gsl::span<T> outputLeft, gsl::span<T> outputRight) noexcept
void readInterleaved(absl::Span<const T> input, absl::Span<T> outputLeft, absl::Span<T> outputRight) noexcept
{
// The size of the output is not big enough for the input...
ASSERT(outputLeft.size() >= input.size() / 2);
@ -20,7 +22,7 @@ void readInterleaved(gsl::span<const T> input, gsl::span<T> outputLeft, gsl::spa
}
template<class T, bool SIMD=SIMDConfig::writeInterleaved>
void writeInterleaved(gsl::span<const T> inputLeft, gsl::span<const T> inputRight, gsl::span<T> output) noexcept
void writeInterleaved(absl::Span<const T> inputLeft, absl::Span<const T> inputRight, absl::Span<T> output) noexcept
{
ASSERT(inputLeft.size() <= output.size() / 2);
ASSERT(inputRight.size() <= output.size() / 2);
@ -37,31 +39,43 @@ void writeInterleaved(gsl::span<const T> inputLeft, gsl::span<const T> inputRigh
// Specializations
template<>
void writeInterleaved<float, true>(gsl::span<const float> inputLeft, gsl::span<const float> inputRight, gsl::span<float> output) noexcept;
void writeInterleaved<float, true>(absl::Span<const float> inputLeft, absl::Span<const float> inputRight, absl::Span<float> output) noexcept;
template<>
void readInterleaved<float, true>(gsl::span<const float> input, gsl::span<float> outputLeft, gsl::span<float> outputRight) noexcept;
void readInterleaved<float, true>(absl::Span<const float> input, absl::Span<float> outputLeft, absl::Span<float> outputRight) noexcept;
template<class T, bool SIMD=SIMDConfig::fill>
void fill(gsl::span<T> output, T value) noexcept
void fill(absl::Span<T> output, T value) noexcept
{
std::fill(output.begin(), output.end(), value);
}
template<>
void fill<float, true>(gsl::span<float> output, float value) noexcept;
void fill<float, true>(absl::Span<float> output, float value) noexcept;
template<class Type, bool SIMD=SIMDConfig::useSIMD>
void exp(absl::Span<const Type> input, absl::Span<Type> output) noexcept
{
ASSERT(output.size() >= input.size());
auto sentinel = std::min(input.size(), output.size());
for (decltype(sentinel) i = 0; i < sentinel; ++i)
output[i] = std::exp(input[i]);
}
template<>
void exp<float, true>(absl::Span<const float> input, absl::Span<float> output) noexcept;
template<class T, bool SIMD=SIMDConfig::useSIMD>
void loopingSFZIndex(gsl::span<const T> inputLeft, gsl::span<const T> inputRight, gsl::span<T> output);
void loopingSFZIndex(absl::Span<const T> inputLeft, absl::Span<const T> inputRight, absl::Span<T> output);
template<class T, bool SIMD=SIMDConfig::useSIMD>
void linearRamp(gsl::span<T> output, T start, T end);
void linearRamp(absl::Span<T> output, T start, T end);
template<class T, bool SIMD=SIMDConfig::useSIMD>
void exponentialRamp(gsl::span<T> output, T start, T end);
void exponentialRamp(absl::Span<T> output, T start, T end);
template<class T, bool SIMD=SIMDConfig::useSIMD>
void applyGain(T gain, gsl::span<T> output);
void applyGain(T gain, absl::Span<T> output);
template<class T, bool SIMD=SIMDConfig::useSIMD>
void applyGain(gsl::span<const T> gain, gsl::span<T> output);
void applyGain(absl::Span<const T> gain, absl::Span<T> output);

View file

@ -1,12 +1,13 @@
#include "SIMDHelpers.h"
#include "Helpers.h"
#include "x86intrin.h"
#include "mathfuns/sse_mathfun.h"
constexpr int TypeAlignment { 4 };
using Type = float;
template<>
void readInterleaved<Type, true>(gsl::span<const Type> input, gsl::span<Type> outputLeft, gsl::span<Type> outputRight) noexcept
void readInterleaved<Type, true>(absl::Span<const Type> input, absl::Span<Type> outputLeft, absl::Span<Type> outputRight) noexcept
{
// The size of the outputs is not big enough for the input...
ASSERT(outputLeft.size() >= input.size() / 2);
@ -36,7 +37,7 @@ void readInterleaved<Type, true>(gsl::span<const Type> input, gsl::span<Type> ou
_mm_storeu_ps(rOut, register1);
lOut += TypeAlignment;
rOut += TypeAlignment;
}
}
inputSentinel = input.end() - 1;
while (in < inputSentinel && lOut < outputLeft.end() && rOut < outputRight.end())
@ -47,7 +48,7 @@ void readInterleaved<Type, true>(gsl::span<const Type> input, gsl::span<Type> ou
}
template<>
void writeInterleaved<Type, true>(gsl::span<const Type> inputLeft, gsl::span<const Type> inputRight, gsl::span<Type> output) noexcept
void writeInterleaved<Type, true>(absl::Span<const Type> inputLeft, absl::Span<const Type> inputRight, absl::Span<Type> output) noexcept
{
// The size of the output is not big enough for the inputs...
ASSERT(inputLeft.size() <= output.size() / 2);
@ -89,7 +90,7 @@ void writeInterleaved<Type, true>(gsl::span<const Type> inputLeft, gsl::span<con
}
template<>
void fill<float, true>(gsl::span<float> output, float value) noexcept
void fill<float, true>(absl::Span<float> output, float value) noexcept
{
const auto mmValue = _mm_set_ps1(value);
auto* out = output.begin();
@ -104,4 +105,19 @@ void fill<float, true>(gsl::span<float> output, float value) noexcept
while (out < output.end())
*out++ = value;
}
template<>
void exp<float, true>(absl::Span<const float> input, absl::Span<float> output) noexcept
{
ASSERT(output.size() >= input.size());
auto* in = input.begin();
auto* out = output.begin();
auto* sentinel = in + std::min(input.size(), output.size());
while (in < sentinel)
{
_mm_storeu_ps(out, exp_ps(_mm_loadu_ps(in)));
out += 4;
in += 4;
}
}

View file

@ -2,7 +2,7 @@
#include "Buffer.h"
#include "Helpers.h"
#include "Globals.h"
#include "gsl/gsl-lite.hpp"
#include <array>
#include <iostream>
#include <type_traits>
@ -50,17 +50,17 @@ public:
void fill(Type value) noexcept
{
::fill<Type>(leftBuffer, value);
::fill<Type>(rightBuffer, value);
::fill<Type>(absl::MakeSpan(leftBuffer), value);
::fill<Type>(absl::MakeSpan(rightBuffer), value);
}
void readInterleaved(gsl::span<const Type> input) noexcept
void readInterleaved(absl::Span<const Type> input) noexcept
{
ASSERT(input.size() <= numChannels * numFrames);
::readInterleaved<Type>(input, leftBuffer, rightBuffer);
::readInterleaved<Type>(input, absl::MakeSpan(leftBuffer), absl::MakeSpan(rightBuffer));
}
void writeInterleaved(gsl::span<Type> output) noexcept
void writeInterleaved(absl::Span<Type> output) noexcept
{
ASSERT(output.size() >= numChannels * numFrames);
::writeInterleaved<Type>(leftBuffer, rightBuffer, output);

View file

@ -0,0 +1,301 @@
/* NEON implementation of sin, cos, exp and log
Inspired by Intel Approximate Math library, and based on the
corresponding algorithms of the cephes math library
*/
/* Copyright (C) 2011 Julien Pommier
This software is provided 'as-is', without any express or implied
warranty. In no event will the authors be held liable for any damages
arising from the use of this software.
Permission is granted to anyone to use this software for any purpose,
including commercial applications, and to alter it and redistribute it
freely, subject to the following restrictions:
1. The origin of this software must not be misrepresented; you must not
claim that you wrote the original software. If you use this software
in a product, an acknowledgment in the product documentation would be
appreciated but is not required.
2. Altered source versions must be plainly marked as such, and must not be
misrepresented as being the original software.
3. This notice may not be removed or altered from any source distribution.
(this is the zlib license)
*/
#include <arm_neon.h>
typedef float32x4_t v4sf; // vector of 4 float
typedef uint32x4_t v4su; // vector of 4 uint32
typedef int32x4_t v4si; // vector of 4 uint32
#define c_inv_mant_mask ~0x7f800000u
#define c_cephes_SQRTHF 0.707106781186547524
#define c_cephes_log_p0 7.0376836292E-2
#define c_cephes_log_p1 - 1.1514610310E-1
#define c_cephes_log_p2 1.1676998740E-1
#define c_cephes_log_p3 - 1.2420140846E-1
#define c_cephes_log_p4 + 1.4249322787E-1
#define c_cephes_log_p5 - 1.6668057665E-1
#define c_cephes_log_p6 + 2.0000714765E-1
#define c_cephes_log_p7 - 2.4999993993E-1
#define c_cephes_log_p8 + 3.3333331174E-1
#define c_cephes_log_q1 -2.12194440e-4
#define c_cephes_log_q2 0.693359375
/* natural logarithm computed for 4 simultaneous float
return NaN for x <= 0
*/
v4sf log_ps(v4sf x) {
v4sf one = vdupq_n_f32(1);
x = vmaxq_f32(x, vdupq_n_f32(0)); /* force flush to zero on denormal values */
v4su invalid_mask = vcleq_f32(x, vdupq_n_f32(0));
v4si ux = vreinterpretq_s32_f32(x);
v4si emm0 = vshrq_n_s32(ux, 23);
/* keep only the fractional part */
ux = vandq_s32(ux, vdupq_n_s32(c_inv_mant_mask));
ux = vorrq_s32(ux, vreinterpretq_s32_f32(vdupq_n_f32(0.5f)));
x = vreinterpretq_f32_s32(ux);
emm0 = vsubq_s32(emm0, vdupq_n_s32(0x7f));
v4sf e = vcvtq_f32_s32(emm0);
e = vaddq_f32(e, one);
/* part2:
if( x < SQRTHF ) {
e -= 1;
x = x + x - 1.0;
} else { x = x - 1.0; }
*/
v4su mask = vcltq_f32(x, vdupq_n_f32(c_cephes_SQRTHF));
v4sf tmp = vreinterpretq_f32_u32(vandq_u32(vreinterpretq_u32_f32(x), mask));
x = vsubq_f32(x, one);
e = vsubq_f32(e, vreinterpretq_f32_u32(vandq_u32(vreinterpretq_u32_f32(one), mask)));
x = vaddq_f32(x, tmp);
v4sf z = vmulq_f32(x,x);
v4sf y = vdupq_n_f32(c_cephes_log_p0);
y = vmulq_f32(y, x);
y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p1));
y = vmulq_f32(y, x);
y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p2));
y = vmulq_f32(y, x);
y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p3));
y = vmulq_f32(y, x);
y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p4));
y = vmulq_f32(y, x);
y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p5));
y = vmulq_f32(y, x);
y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p6));
y = vmulq_f32(y, x);
y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p7));
y = vmulq_f32(y, x);
y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p8));
y = vmulq_f32(y, x);
y = vmulq_f32(y, z);
tmp = vmulq_f32(e, vdupq_n_f32(c_cephes_log_q1));
y = vaddq_f32(y, tmp);
tmp = vmulq_f32(z, vdupq_n_f32(0.5f));
y = vsubq_f32(y, tmp);
tmp = vmulq_f32(e, vdupq_n_f32(c_cephes_log_q2));
x = vaddq_f32(x, y);
x = vaddq_f32(x, tmp);
x = vreinterpretq_f32_u32(vorrq_u32(vreinterpretq_u32_f32(x), invalid_mask)); // negative arg will be NAN
return x;
}
#define c_exp_hi 88.3762626647949f
#define c_exp_lo -88.3762626647949f
#define c_cephes_LOG2EF 1.44269504088896341
#define c_cephes_exp_C1 0.693359375
#define c_cephes_exp_C2 -2.12194440e-4
#define c_cephes_exp_p0 1.9875691500E-4
#define c_cephes_exp_p1 1.3981999507E-3
#define c_cephes_exp_p2 8.3334519073E-3
#define c_cephes_exp_p3 4.1665795894E-2
#define c_cephes_exp_p4 1.6666665459E-1
#define c_cephes_exp_p5 5.0000001201E-1
/* exp() computed for 4 float at once */
v4sf exp_ps(v4sf x) {
v4sf tmp, fx;
v4sf one = vdupq_n_f32(1);
x = vminq_f32(x, vdupq_n_f32(c_exp_hi));
x = vmaxq_f32(x, vdupq_n_f32(c_exp_lo));
/* express exp(x) as exp(g + n*log(2)) */
fx = vmlaq_f32(vdupq_n_f32(0.5f), x, vdupq_n_f32(c_cephes_LOG2EF));
/* perform a floorf */
tmp = vcvtq_f32_s32(vcvtq_s32_f32(fx));
/* if greater, substract 1 */
v4su mask = vcgtq_f32(tmp, fx);
mask = vandq_u32(mask, vreinterpretq_u32_f32(one));
fx = vsubq_f32(tmp, vreinterpretq_f32_u32(mask));
tmp = vmulq_f32(fx, vdupq_n_f32(c_cephes_exp_C1));
v4sf z = vmulq_f32(fx, vdupq_n_f32(c_cephes_exp_C2));
x = vsubq_f32(x, tmp);
x = vsubq_f32(x, z);
static const float cephes_exp_p[6] = { c_cephes_exp_p0, c_cephes_exp_p1, c_cephes_exp_p2, c_cephes_exp_p3, c_cephes_exp_p4, c_cephes_exp_p5 };
v4sf y = vld1q_dup_f32(cephes_exp_p+0);
v4sf c1 = vld1q_dup_f32(cephes_exp_p+1);
v4sf c2 = vld1q_dup_f32(cephes_exp_p+2);
v4sf c3 = vld1q_dup_f32(cephes_exp_p+3);
v4sf c4 = vld1q_dup_f32(cephes_exp_p+4);
v4sf c5 = vld1q_dup_f32(cephes_exp_p+5);
y = vmulq_f32(y, x);
z = vmulq_f32(x,x);
y = vaddq_f32(y, c1);
y = vmulq_f32(y, x);
y = vaddq_f32(y, c2);
y = vmulq_f32(y, x);
y = vaddq_f32(y, c3);
y = vmulq_f32(y, x);
y = vaddq_f32(y, c4);
y = vmulq_f32(y, x);
y = vaddq_f32(y, c5);
y = vmulq_f32(y, z);
y = vaddq_f32(y, x);
y = vaddq_f32(y, one);
/* build 2^n */
int32x4_t mm;
mm = vcvtq_s32_f32(fx);
mm = vaddq_s32(mm, vdupq_n_s32(0x7f));
mm = vshlq_n_s32(mm, 23);
v4sf pow2n = vreinterpretq_f32_s32(mm);
y = vmulq_f32(y, pow2n);
return y;
}
#define c_minus_cephes_DP1 -0.78515625
#define c_minus_cephes_DP2 -2.4187564849853515625e-4
#define c_minus_cephes_DP3 -3.77489497744594108e-8
#define c_sincof_p0 -1.9515295891E-4
#define c_sincof_p1 8.3321608736E-3
#define c_sincof_p2 -1.6666654611E-1
#define c_coscof_p0 2.443315711809948E-005
#define c_coscof_p1 -1.388731625493765E-003
#define c_coscof_p2 4.166664568298827E-002
#define c_cephes_FOPI 1.27323954473516 // 4 / M_PI
/* evaluation of 4 sines & cosines at once.
The code is the exact rewriting of the cephes sinf function.
Precision is excellent as long as x < 8192 (I did not bother to
take into account the special handling they have for greater values
-- it does not return garbage for arguments over 8192, though, but
the extra precision is missing).
Note that it is such that sinf((float)M_PI) = 8.74e-8, which is the
surprising but correct result.
Note also that when you compute sin(x), cos(x) is available at
almost no extra price so both sin_ps and cos_ps make use of
sincos_ps..
*/
void sincos_ps(v4sf x, v4sf *ysin, v4sf *ycos) { // any x
v4sf xmm1, xmm2, xmm3, y;
v4su emm2;
v4su sign_mask_sin, sign_mask_cos;
sign_mask_sin = vcltq_f32(x, vdupq_n_f32(0));
x = vabsq_f32(x);
/* scale by 4/Pi */
y = vmulq_f32(x, vdupq_n_f32(c_cephes_FOPI));
/* store the integer part of y in mm0 */
emm2 = vcvtq_u32_f32(y);
/* j=(j+1) & (~1) (see the cephes sources) */
emm2 = vaddq_u32(emm2, vdupq_n_u32(1));
emm2 = vandq_u32(emm2, vdupq_n_u32(~1));
y = vcvtq_f32_u32(emm2);
/* get the polynom selection mask
there is one polynom for 0 <= x <= Pi/4
and another one for Pi/4<x<=Pi/2
Both branches will be computed.
*/
v4su poly_mask = vtstq_u32(emm2, vdupq_n_u32(2));
/* The magic pass: "Extended precision modular arithmetic"
x = ((x - y * DP1) - y * DP2) - y * DP3; */
xmm1 = vmulq_n_f32(y, c_minus_cephes_DP1);
xmm2 = vmulq_n_f32(y, c_minus_cephes_DP2);
xmm3 = vmulq_n_f32(y, c_minus_cephes_DP3);
x = vaddq_f32(x, xmm1);
x = vaddq_f32(x, xmm2);
x = vaddq_f32(x, xmm3);
sign_mask_sin = veorq_u32(sign_mask_sin, vtstq_u32(emm2, vdupq_n_u32(4)));
sign_mask_cos = vtstq_u32(vsubq_u32(emm2, vdupq_n_u32(2)), vdupq_n_u32(4));
/* Evaluate the first polynom (0 <= x <= Pi/4) in y1,
and the second polynom (Pi/4 <= x <= 0) in y2 */
v4sf z = vmulq_f32(x,x);
v4sf y1, y2;
y1 = vmulq_n_f32(z, c_coscof_p0);
y2 = vmulq_n_f32(z, c_sincof_p0);
y1 = vaddq_f32(y1, vdupq_n_f32(c_coscof_p1));
y2 = vaddq_f32(y2, vdupq_n_f32(c_sincof_p1));
y1 = vmulq_f32(y1, z);
y2 = vmulq_f32(y2, z);
y1 = vaddq_f32(y1, vdupq_n_f32(c_coscof_p2));
y2 = vaddq_f32(y2, vdupq_n_f32(c_sincof_p2));
y1 = vmulq_f32(y1, z);
y2 = vmulq_f32(y2, z);
y1 = vmulq_f32(y1, z);
y2 = vmulq_f32(y2, x);
y1 = vsubq_f32(y1, vmulq_f32(z, vdupq_n_f32(0.5f)));
y2 = vaddq_f32(y2, x);
y1 = vaddq_f32(y1, vdupq_n_f32(1));
/* select the correct result from the two polynoms */
v4sf ys = vbslq_f32(poly_mask, y1, y2);
v4sf yc = vbslq_f32(poly_mask, y2, y1);
*ysin = vbslq_f32(sign_mask_sin, vnegq_f32(ys), ys);
*ycos = vbslq_f32(sign_mask_cos, yc, vnegq_f32(yc));
}
v4sf sin_ps(v4sf x) {
v4sf ysin, ycos;
sincos_ps(x, &ysin, &ycos);
return ysin;
}
v4sf cos_ps(v4sf x) {
v4sf ysin, ycos;
sincos_ps(x, &ysin, &ycos);
return ycos;
}

View file

@ -0,0 +1,713 @@
/* SIMD (SSE1+MMX or SSE2) implementation of sin, cos, exp and log
Inspired by Intel Approximate Math library, and based on the
corresponding algorithms of the cephes math library
The default is to use the SSE1 version. If you define USE_SSE2 the
the SSE2 intrinsics will be used in place of the MMX intrinsics. Do
not expect any significant performance improvement with SSE2.
*/
/* Copyright (C) 2007 Julien Pommier
This software is provided 'as-is', without any express or implied
warranty. In no event will the authors be held liable for any damages
arising from the use of this software.
Permission is granted to anyone to use this software for any purpose,
including commercial applications, and to alter it and redistribute it
freely, subject to the following restrictions:
1. The origin of this software must not be misrepresented; you must not
claim that you wrote the original software. If you use this software
in a product, an acknowledgment in the product documentation would be
appreciated but is not required.
2. Altered source versions must be plainly marked as such, and must not be
misrepresented as being the original software.
3. This notice may not be removed or altered from any source distribution.
(this is the zlib license)
*/
#include <xmmintrin.h>
/* yes I know, the top of this file is quite ugly */
#ifdef _MSC_VER /* visual c++ */
# define ALIGN16_BEG __declspec(align(16))
# define ALIGN16_END
#else /* gcc or icc */
# define ALIGN16_BEG
# define ALIGN16_END __attribute__((aligned(16)))
#endif
#define USE_SSE2
/* __m128 is ugly to write */
typedef __m128 v4sf; // vector of 4 float (sse1)
#ifdef USE_SSE2
# include <emmintrin.h>
typedef __m128i v4si; // vector of 4 int (sse2)
#else
typedef __m64 v2si; // vector of 2 int (mmx)
#endif
/* declare some SSE constants -- why can't I figure a better way to do that? */
#define _PS_CONST(Name, Val) \
static const ALIGN16_BEG float _ps_##Name[4] ALIGN16_END = { Val, Val, Val, Val }
#define _PI32_CONST(Name, Val) \
static const ALIGN16_BEG int _pi32_##Name[4] ALIGN16_END = { Val, Val, Val, Val }
#define _PS_CONST_TYPE(Name, Type, Val) \
static const ALIGN16_BEG Type _ps_##Name[4] ALIGN16_END = { Val, Val, Val, Val }
_PS_CONST(1 , 1.0f);
_PS_CONST(0p5, 0.5f);
/* the smallest non denormalized float number */
_PS_CONST_TYPE(min_norm_pos, int, 0x00800000);
_PS_CONST_TYPE(mant_mask, int, 0x7f800000);
_PS_CONST_TYPE(inv_mant_mask, int, ~0x7f800000);
_PS_CONST_TYPE(sign_mask, int, (int)0x80000000);
_PS_CONST_TYPE(inv_sign_mask, int, ~0x80000000);
_PI32_CONST(1, 1);
_PI32_CONST(inv1, ~1);
_PI32_CONST(2, 2);
_PI32_CONST(4, 4);
_PI32_CONST(0x7f, 0x7f);
_PS_CONST(cephes_SQRTHF, 0.707106781186547524);
_PS_CONST(cephes_log_p0, 7.0376836292E-2);
_PS_CONST(cephes_log_p1, - 1.1514610310E-1);
_PS_CONST(cephes_log_p2, 1.1676998740E-1);
_PS_CONST(cephes_log_p3, - 1.2420140846E-1);
_PS_CONST(cephes_log_p4, + 1.4249322787E-1);
_PS_CONST(cephes_log_p5, - 1.6668057665E-1);
_PS_CONST(cephes_log_p6, + 2.0000714765E-1);
_PS_CONST(cephes_log_p7, - 2.4999993993E-1);
_PS_CONST(cephes_log_p8, + 3.3333331174E-1);
_PS_CONST(cephes_log_q1, -2.12194440e-4);
_PS_CONST(cephes_log_q2, 0.693359375);
#ifndef USE_SSE2
typedef union xmm_mm_union {
__m128 xmm;
__m64 mm[2];
} xmm_mm_union;
#define COPY_XMM_TO_MM(xmm_, mm0_, mm1_) { \
xmm_mm_union u; u.xmm = xmm_; \
mm0_ = u.mm[0]; \
mm1_ = u.mm[1]; \
}
#define COPY_MM_TO_XMM(mm0_, mm1_, xmm_) { \
xmm_mm_union u; u.mm[0]=mm0_; u.mm[1]=mm1_; xmm_ = u.xmm; \
}
#endif // USE_SSE2
/* natural logarithm computed for 4 simultaneous float
return NaN for x <= 0
*/
v4sf log_ps(v4sf x) {
#ifdef USE_SSE2
v4si emm0;
#else
v2si mm0, mm1;
#endif
v4sf one = *(v4sf*)_ps_1;
v4sf invalid_mask = _mm_cmple_ps(x, _mm_setzero_ps());
x = _mm_max_ps(x, *(v4sf*)_ps_min_norm_pos); /* cut off denormalized stuff */
#ifndef USE_SSE2
/* part 1: x = frexpf(x, &e); */
COPY_XMM_TO_MM(x, mm0, mm1);
mm0 = _mm_srli_pi32(mm0, 23);
mm1 = _mm_srli_pi32(mm1, 23);
#else
emm0 = _mm_srli_epi32(_mm_castps_si128(x), 23);
#endif
/* keep only the fractional part */
x = _mm_and_ps(x, *(v4sf*)_ps_inv_mant_mask);
x = _mm_or_ps(x, *(v4sf*)_ps_0p5);
#ifndef USE_SSE2
/* now e=mm0:mm1 contain the really base-2 exponent */
mm0 = _mm_sub_pi32(mm0, *(v2si*)_pi32_0x7f);
mm1 = _mm_sub_pi32(mm1, *(v2si*)_pi32_0x7f);
v4sf e = _mm_cvtpi32x2_ps(mm0, mm1);
_mm_empty(); /* bye bye mmx */
#else
emm0 = _mm_sub_epi32(emm0, *(v4si*)_pi32_0x7f);
v4sf e = _mm_cvtepi32_ps(emm0);
#endif
e = _mm_add_ps(e, one);
/* part2:
if( x < SQRTHF ) {
e -= 1;
x = x + x - 1.0;
} else { x = x - 1.0; }
*/
v4sf mask = _mm_cmplt_ps(x, *(v4sf*)_ps_cephes_SQRTHF);
v4sf tmp = _mm_and_ps(x, mask);
x = _mm_sub_ps(x, one);
e = _mm_sub_ps(e, _mm_and_ps(one, mask));
x = _mm_add_ps(x, tmp);
v4sf z = _mm_mul_ps(x,x);
v4sf y = *(v4sf*)_ps_cephes_log_p0;
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p1);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p2);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p3);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p4);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p5);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p6);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p7);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p8);
y = _mm_mul_ps(y, x);
y = _mm_mul_ps(y, z);
tmp = _mm_mul_ps(e, *(v4sf*)_ps_cephes_log_q1);
y = _mm_add_ps(y, tmp);
tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
y = _mm_sub_ps(y, tmp);
tmp = _mm_mul_ps(e, *(v4sf*)_ps_cephes_log_q2);
x = _mm_add_ps(x, y);
x = _mm_add_ps(x, tmp);
x = _mm_or_ps(x, invalid_mask); // negative arg will be NAN
return x;
}
_PS_CONST(exp_hi, 88.3762626647949f);
_PS_CONST(exp_lo, -88.3762626647949f);
_PS_CONST(cephes_LOG2EF, 1.44269504088896341);
_PS_CONST(cephes_exp_C1, 0.693359375);
_PS_CONST(cephes_exp_C2, -2.12194440e-4);
_PS_CONST(cephes_exp_p0, 1.9875691500E-4);
_PS_CONST(cephes_exp_p1, 1.3981999507E-3);
_PS_CONST(cephes_exp_p2, 8.3334519073E-3);
_PS_CONST(cephes_exp_p3, 4.1665795894E-2);
_PS_CONST(cephes_exp_p4, 1.6666665459E-1);
_PS_CONST(cephes_exp_p5, 5.0000001201E-1);
v4sf exp_ps(v4sf x) {
v4sf tmp = _mm_setzero_ps(), fx;
#ifdef USE_SSE2
v4si emm0;
#else
v2si mm0, mm1;
#endif
v4sf one = *(v4sf*)_ps_1;
x = _mm_min_ps(x, *(v4sf*)_ps_exp_hi);
x = _mm_max_ps(x, *(v4sf*)_ps_exp_lo);
/* express exp(x) as exp(g + n*log(2)) */
fx = _mm_mul_ps(x, *(v4sf*)_ps_cephes_LOG2EF);
fx = _mm_add_ps(fx, *(v4sf*)_ps_0p5);
/* how to perform a floorf with SSE: just below */
#ifndef USE_SSE2
/* step 1 : cast to int */
tmp = _mm_movehl_ps(tmp, fx);
mm0 = _mm_cvttps_pi32(fx);
mm1 = _mm_cvttps_pi32(tmp);
/* step 2 : cast back to float */
tmp = _mm_cvtpi32x2_ps(mm0, mm1);
#else
emm0 = _mm_cvttps_epi32(fx);
tmp = _mm_cvtepi32_ps(emm0);
#endif
/* if greater, substract 1 */
v4sf mask = _mm_cmpgt_ps(tmp, fx);
mask = _mm_and_ps(mask, one);
fx = _mm_sub_ps(tmp, mask);
tmp = _mm_mul_ps(fx, *(v4sf*)_ps_cephes_exp_C1);
v4sf z = _mm_mul_ps(fx, *(v4sf*)_ps_cephes_exp_C2);
x = _mm_sub_ps(x, tmp);
x = _mm_sub_ps(x, z);
z = _mm_mul_ps(x,x);
v4sf y = *(v4sf*)_ps_cephes_exp_p0;
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p1);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p2);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p3);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p4);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p5);
y = _mm_mul_ps(y, z);
y = _mm_add_ps(y, x);
y = _mm_add_ps(y, one);
/* build 2^n */
#ifndef USE_SSE2
z = _mm_movehl_ps(z, fx);
mm0 = _mm_cvttps_pi32(fx);
mm1 = _mm_cvttps_pi32(z);
mm0 = _mm_add_pi32(mm0, *(v2si*)_pi32_0x7f);
mm1 = _mm_add_pi32(mm1, *(v2si*)_pi32_0x7f);
mm0 = _mm_slli_pi32(mm0, 23);
mm1 = _mm_slli_pi32(mm1, 23);
v4sf pow2n;
COPY_MM_TO_XMM(mm0, mm1, pow2n);
_mm_empty();
#else
emm0 = _mm_cvttps_epi32(fx);
emm0 = _mm_add_epi32(emm0, *(v4si*)_pi32_0x7f);
emm0 = _mm_slli_epi32(emm0, 23);
v4sf pow2n = _mm_castsi128_ps(emm0);
#endif
y = _mm_mul_ps(y, pow2n);
return y;
}
_PS_CONST(minus_cephes_DP1, -0.78515625);
_PS_CONST(minus_cephes_DP2, -2.4187564849853515625e-4);
_PS_CONST(minus_cephes_DP3, -3.77489497744594108e-8);
_PS_CONST(sincof_p0, -1.9515295891E-4);
_PS_CONST(sincof_p1, 8.3321608736E-3);
_PS_CONST(sincof_p2, -1.6666654611E-1);
_PS_CONST(coscof_p0, 2.443315711809948E-005);
_PS_CONST(coscof_p1, -1.388731625493765E-003);
_PS_CONST(coscof_p2, 4.166664568298827E-002);
_PS_CONST(cephes_FOPI, 1.27323954473516); // 4 / M_PI
/* evaluation of 4 sines at onces, using only SSE1+MMX intrinsics so
it runs also on old athlons XPs and the pentium III of your grand
mother.
The code is the exact rewriting of the cephes sinf function.
Precision is excellent as long as x < 8192 (I did not bother to
take into account the special handling they have for greater values
-- it does not return garbage for arguments over 8192, though, but
the extra precision is missing).
Note that it is such that sinf((float)M_PI) = 8.74e-8, which is the
surprising but correct result.
Performance is also surprisingly good, 1.33 times faster than the
macos vsinf SSE2 function, and 1.5 times faster than the
__vrs4_sinf of amd's ACML (which is only available in 64 bits). Not
too bad for an SSE1 function (with no special tuning) !
However the latter libraries probably have a much better handling of NaN,
Inf, denormalized and other special arguments..
On my core 1 duo, the execution of this function takes approximately 95 cycles.
From what I have observed on the experiments with Intel AMath lib, switching to an
SSE2 version would improve the perf by only 10%.
Since it is based on SSE intrinsics, it has to be compiled at -O2 to
deliver full speed.
*/
v4sf sin_ps(v4sf x) { // any x
v4sf xmm1, xmm2 = _mm_setzero_ps(), xmm3, sign_bit, y;
#ifdef USE_SSE2
v4si emm0, emm2;
#else
v2si mm0, mm1, mm2, mm3;
#endif
sign_bit = x;
/* take the absolute value */
x = _mm_and_ps(x, *(v4sf*)_ps_inv_sign_mask);
/* extract the sign bit (upper one) */
sign_bit = _mm_and_ps(sign_bit, *(v4sf*)_ps_sign_mask);
/* scale by 4/Pi */
y = _mm_mul_ps(x, *(v4sf*)_ps_cephes_FOPI);
#ifdef USE_SSE2
/* store the integer part of y in mm0 */
emm2 = _mm_cvttps_epi32(y);
/* j=(j+1) & (~1) (see the cephes sources) */
emm2 = _mm_add_epi32(emm2, *(v4si*)_pi32_1);
emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_inv1);
y = _mm_cvtepi32_ps(emm2);
/* get the swap sign flag */
emm0 = _mm_and_si128(emm2, *(v4si*)_pi32_4);
emm0 = _mm_slli_epi32(emm0, 29);
/* get the polynom selection mask
there is one polynom for 0 <= x <= Pi/4
and another one for Pi/4<x<=Pi/2
Both branches will be computed.
*/
emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_2);
emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
v4sf swap_sign_bit = _mm_castsi128_ps(emm0);
v4sf poly_mask = _mm_castsi128_ps(emm2);
sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);
#else
/* store the integer part of y in mm0:mm1 */
xmm2 = _mm_movehl_ps(xmm2, y);
mm2 = _mm_cvttps_pi32(y);
mm3 = _mm_cvttps_pi32(xmm2);
/* j=(j+1) & (~1) (see the cephes sources) */
mm2 = _mm_add_pi32(mm2, *(v2si*)_pi32_1);
mm3 = _mm_add_pi32(mm3, *(v2si*)_pi32_1);
mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_inv1);
mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_inv1);
y = _mm_cvtpi32x2_ps(mm2, mm3);
/* get the swap sign flag */
mm0 = _mm_and_si64(mm2, *(v2si*)_pi32_4);
mm1 = _mm_and_si64(mm3, *(v2si*)_pi32_4);
mm0 = _mm_slli_pi32(mm0, 29);
mm1 = _mm_slli_pi32(mm1, 29);
/* get the polynom selection mask */
mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_2);
mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_2);
mm2 = _mm_cmpeq_pi32(mm2, _mm_setzero_si64());
mm3 = _mm_cmpeq_pi32(mm3, _mm_setzero_si64());
v4sf swap_sign_bit, poly_mask;
COPY_MM_TO_XMM(mm0, mm1, swap_sign_bit);
COPY_MM_TO_XMM(mm2, mm3, poly_mask);
sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);
_mm_empty(); /* good-bye mmx */
#endif
/* The magic pass: "Extended precision modular arithmetic"
x = ((x - y * DP1) - y * DP2) - y * DP3; */
xmm1 = *(v4sf*)_ps_minus_cephes_DP1;
xmm2 = *(v4sf*)_ps_minus_cephes_DP2;
xmm3 = *(v4sf*)_ps_minus_cephes_DP3;
xmm1 = _mm_mul_ps(y, xmm1);
xmm2 = _mm_mul_ps(y, xmm2);
xmm3 = _mm_mul_ps(y, xmm3);
x = _mm_add_ps(x, xmm1);
x = _mm_add_ps(x, xmm2);
x = _mm_add_ps(x, xmm3);
/* Evaluate the first polynom (0 <= x <= Pi/4) */
y = *(v4sf*)_ps_coscof_p0;
v4sf z = _mm_mul_ps(x,x);
y = _mm_mul_ps(y, z);
y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p1);
y = _mm_mul_ps(y, z);
y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p2);
y = _mm_mul_ps(y, z);
y = _mm_mul_ps(y, z);
v4sf tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
y = _mm_sub_ps(y, tmp);
y = _mm_add_ps(y, *(v4sf*)_ps_1);
/* Evaluate the second polynom (Pi/4 <= x <= 0) */
v4sf y2 = *(v4sf*)_ps_sincof_p0;
y2 = _mm_mul_ps(y2, z);
y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p1);
y2 = _mm_mul_ps(y2, z);
y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p2);
y2 = _mm_mul_ps(y2, z);
y2 = _mm_mul_ps(y2, x);
y2 = _mm_add_ps(y2, x);
/* select the correct result from the two polynoms */
xmm3 = poly_mask;
y2 = _mm_and_ps(xmm3, y2); //, xmm3);
y = _mm_andnot_ps(xmm3, y);
y = _mm_add_ps(y,y2);
/* update the sign */
y = _mm_xor_ps(y, sign_bit);
return y;
}
/* almost the same as sin_ps */
v4sf cos_ps(v4sf x) { // any x
v4sf xmm1, xmm2 = _mm_setzero_ps(), xmm3, y;
#ifdef USE_SSE2
v4si emm0, emm2;
#else
v2si mm0, mm1, mm2, mm3;
#endif
/* take the absolute value */
x = _mm_and_ps(x, *(v4sf*)_ps_inv_sign_mask);
/* scale by 4/Pi */
y = _mm_mul_ps(x, *(v4sf*)_ps_cephes_FOPI);
#ifdef USE_SSE2
/* store the integer part of y in mm0 */
emm2 = _mm_cvttps_epi32(y);
/* j=(j+1) & (~1) (see the cephes sources) */
emm2 = _mm_add_epi32(emm2, *(v4si*)_pi32_1);
emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_inv1);
y = _mm_cvtepi32_ps(emm2);
emm2 = _mm_sub_epi32(emm2, *(v4si*)_pi32_2);
/* get the swap sign flag */
emm0 = _mm_andnot_si128(emm2, *(v4si*)_pi32_4);
emm0 = _mm_slli_epi32(emm0, 29);
/* get the polynom selection mask */
emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_2);
emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
v4sf sign_bit = _mm_castsi128_ps(emm0);
v4sf poly_mask = _mm_castsi128_ps(emm2);
#else
/* store the integer part of y in mm0:mm1 */
xmm2 = _mm_movehl_ps(xmm2, y);
mm2 = _mm_cvttps_pi32(y);
mm3 = _mm_cvttps_pi32(xmm2);
/* j=(j+1) & (~1) (see the cephes sources) */
mm2 = _mm_add_pi32(mm2, *(v2si*)_pi32_1);
mm3 = _mm_add_pi32(mm3, *(v2si*)_pi32_1);
mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_inv1);
mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_inv1);
y = _mm_cvtpi32x2_ps(mm2, mm3);
mm2 = _mm_sub_pi32(mm2, *(v2si*)_pi32_2);
mm3 = _mm_sub_pi32(mm3, *(v2si*)_pi32_2);
/* get the swap sign flag in mm0:mm1 and the
polynom selection mask in mm2:mm3 */
mm0 = _mm_andnot_si64(mm2, *(v2si*)_pi32_4);
mm1 = _mm_andnot_si64(mm3, *(v2si*)_pi32_4);
mm0 = _mm_slli_pi32(mm0, 29);
mm1 = _mm_slli_pi32(mm1, 29);
mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_2);
mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_2);
mm2 = _mm_cmpeq_pi32(mm2, _mm_setzero_si64());
mm3 = _mm_cmpeq_pi32(mm3, _mm_setzero_si64());
v4sf sign_bit, poly_mask;
COPY_MM_TO_XMM(mm0, mm1, sign_bit);
COPY_MM_TO_XMM(mm2, mm3, poly_mask);
_mm_empty(); /* good-bye mmx */
#endif
/* The magic pass: "Extended precision modular arithmetic"
x = ((x - y * DP1) - y * DP2) - y * DP3; */
xmm1 = *(v4sf*)_ps_minus_cephes_DP1;
xmm2 = *(v4sf*)_ps_minus_cephes_DP2;
xmm3 = *(v4sf*)_ps_minus_cephes_DP3;
xmm1 = _mm_mul_ps(y, xmm1);
xmm2 = _mm_mul_ps(y, xmm2);
xmm3 = _mm_mul_ps(y, xmm3);
x = _mm_add_ps(x, xmm1);
x = _mm_add_ps(x, xmm2);
x = _mm_add_ps(x, xmm3);
/* Evaluate the first polynom (0 <= x <= Pi/4) */
y = *(v4sf*)_ps_coscof_p0;
v4sf z = _mm_mul_ps(x,x);
y = _mm_mul_ps(y, z);
y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p1);
y = _mm_mul_ps(y, z);
y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p2);
y = _mm_mul_ps(y, z);
y = _mm_mul_ps(y, z);
v4sf tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
y = _mm_sub_ps(y, tmp);
y = _mm_add_ps(y, *(v4sf*)_ps_1);
/* Evaluate the second polynom (Pi/4 <= x <= 0) */
v4sf y2 = *(v4sf*)_ps_sincof_p0;
y2 = _mm_mul_ps(y2, z);
y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p1);
y2 = _mm_mul_ps(y2, z);
y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p2);
y2 = _mm_mul_ps(y2, z);
y2 = _mm_mul_ps(y2, x);
y2 = _mm_add_ps(y2, x);
/* select the correct result from the two polynoms */
xmm3 = poly_mask;
y2 = _mm_and_ps(xmm3, y2); //, xmm3);
y = _mm_andnot_ps(xmm3, y);
y = _mm_add_ps(y,y2);
/* update the sign */
y = _mm_xor_ps(y, sign_bit);
return y;
}
/* since sin_ps and cos_ps are almost identical, sincos_ps could replace both of them..
it is almost as fast, and gives you a free cosine with your sine */
void sincos_ps(v4sf x, v4sf *s, v4sf *c) {
v4sf xmm1, xmm2, xmm3 = _mm_setzero_ps(), sign_bit_sin, y;
#ifdef USE_SSE2
v4si emm0, emm2, emm4;
#else
v2si mm0, mm1, mm2, mm3, mm4, mm5;
#endif
sign_bit_sin = x;
/* take the absolute value */
x = _mm_and_ps(x, *(v4sf*)_ps_inv_sign_mask);
/* extract the sign bit (upper one) */
sign_bit_sin = _mm_and_ps(sign_bit_sin, *(v4sf*)_ps_sign_mask);
/* scale by 4/Pi */
y = _mm_mul_ps(x, *(v4sf*)_ps_cephes_FOPI);
#ifdef USE_SSE2
/* store the integer part of y in emm2 */
emm2 = _mm_cvttps_epi32(y);
/* j=(j+1) & (~1) (see the cephes sources) */
emm2 = _mm_add_epi32(emm2, *(v4si*)_pi32_1);
emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_inv1);
y = _mm_cvtepi32_ps(emm2);
emm4 = emm2;
/* get the swap sign flag for the sine */
emm0 = _mm_and_si128(emm2, *(v4si*)_pi32_4);
emm0 = _mm_slli_epi32(emm0, 29);
v4sf swap_sign_bit_sin = _mm_castsi128_ps(emm0);
/* get the polynom selection mask for the sine*/
emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_2);
emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
v4sf poly_mask = _mm_castsi128_ps(emm2);
#else
/* store the integer part of y in mm2:mm3 */
xmm3 = _mm_movehl_ps(xmm3, y);
mm2 = _mm_cvttps_pi32(y);
mm3 = _mm_cvttps_pi32(xmm3);
/* j=(j+1) & (~1) (see the cephes sources) */
mm2 = _mm_add_pi32(mm2, *(v2si*)_pi32_1);
mm3 = _mm_add_pi32(mm3, *(v2si*)_pi32_1);
mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_inv1);
mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_inv1);
y = _mm_cvtpi32x2_ps(mm2, mm3);
mm4 = mm2;
mm5 = mm3;
/* get the swap sign flag for the sine */
mm0 = _mm_and_si64(mm2, *(v2si*)_pi32_4);
mm1 = _mm_and_si64(mm3, *(v2si*)_pi32_4);
mm0 = _mm_slli_pi32(mm0, 29);
mm1 = _mm_slli_pi32(mm1, 29);
v4sf swap_sign_bit_sin;
COPY_MM_TO_XMM(mm0, mm1, swap_sign_bit_sin);
/* get the polynom selection mask for the sine */
mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_2);
mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_2);
mm2 = _mm_cmpeq_pi32(mm2, _mm_setzero_si64());
mm3 = _mm_cmpeq_pi32(mm3, _mm_setzero_si64());
v4sf poly_mask;
COPY_MM_TO_XMM(mm2, mm3, poly_mask);
#endif
/* The magic pass: "Extended precision modular arithmetic"
x = ((x - y * DP1) - y * DP2) - y * DP3; */
xmm1 = *(v4sf*)_ps_minus_cephes_DP1;
xmm2 = *(v4sf*)_ps_minus_cephes_DP2;
xmm3 = *(v4sf*)_ps_minus_cephes_DP3;
xmm1 = _mm_mul_ps(y, xmm1);
xmm2 = _mm_mul_ps(y, xmm2);
xmm3 = _mm_mul_ps(y, xmm3);
x = _mm_add_ps(x, xmm1);
x = _mm_add_ps(x, xmm2);
x = _mm_add_ps(x, xmm3);
#ifdef USE_SSE2
emm4 = _mm_sub_epi32(emm4, *(v4si*)_pi32_2);
emm4 = _mm_andnot_si128(emm4, *(v4si*)_pi32_4);
emm4 = _mm_slli_epi32(emm4, 29);
v4sf sign_bit_cos = _mm_castsi128_ps(emm4);
#else
/* get the sign flag for the cosine */
mm4 = _mm_sub_pi32(mm4, *(v2si*)_pi32_2);
mm5 = _mm_sub_pi32(mm5, *(v2si*)_pi32_2);
mm4 = _mm_andnot_si64(mm4, *(v2si*)_pi32_4);
mm5 = _mm_andnot_si64(mm5, *(v2si*)_pi32_4);
mm4 = _mm_slli_pi32(mm4, 29);
mm5 = _mm_slli_pi32(mm5, 29);
v4sf sign_bit_cos;
COPY_MM_TO_XMM(mm4, mm5, sign_bit_cos);
_mm_empty(); /* good-bye mmx */
#endif
sign_bit_sin = _mm_xor_ps(sign_bit_sin, swap_sign_bit_sin);
/* Evaluate the first polynom (0 <= x <= Pi/4) */
v4sf z = _mm_mul_ps(x,x);
y = *(v4sf*)_ps_coscof_p0;
y = _mm_mul_ps(y, z);
y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p1);
y = _mm_mul_ps(y, z);
y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p2);
y = _mm_mul_ps(y, z);
y = _mm_mul_ps(y, z);
v4sf tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
y = _mm_sub_ps(y, tmp);
y = _mm_add_ps(y, *(v4sf*)_ps_1);
/* Evaluate the second polynom (Pi/4 <= x <= 0) */
v4sf y2 = *(v4sf*)_ps_sincof_p0;
y2 = _mm_mul_ps(y2, z);
y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p1);
y2 = _mm_mul_ps(y2, z);
y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p2);
y2 = _mm_mul_ps(y2, z);
y2 = _mm_mul_ps(y2, x);
y2 = _mm_add_ps(y2, x);
/* select the correct result from the two polynoms */
xmm3 = poly_mask;
v4sf ysin2 = _mm_and_ps(xmm3, y2);
v4sf ysin1 = _mm_andnot_ps(xmm3, y);
y2 = _mm_sub_ps(y2,ysin2);
y = _mm_sub_ps(y, ysin1);
xmm1 = _mm_add_ps(ysin1,ysin2);
xmm2 = _mm_add_ps(y,y2);
/* update the sign */
*s = _mm_xor_ps(xmm1, sign_bit_sin);
*c = _mm_xor_ps(xmm2, sign_bit_cos);
}

View file

@ -1,10 +1,10 @@
#include "../sources/OnePoleFilter.h"
#include "catch2/catch.hpp"
#include "cnpy.h"
#include "gsl/gsl-lite.hpp"
#include <string>
#include <filesystem>
#include <algorithm>
#include <absl/types/span.h>
using namespace Catch::literals;
template<class Type>
@ -28,11 +28,11 @@ void testLowpass(const std::filesystem::path& inputNumpyFile, const std::filesys
{
const auto input = cnpy::npy_load(inputNumpyFile.string());
REQUIRE( input.word_size == 8 );
const auto inputSpan = gsl::make_span(input.data<double>(), input.shape[0]);
const auto inputSpan = absl::MakeSpan(input.data<double>(), input.shape[0]);
const auto output = cnpy::npy_load(outputNumpyFile.string());
REQUIRE( output.word_size == 8 );
const auto outputSpan = gsl::make_span(output.data<double>(), output.shape[0]);
const auto outputSpan = absl::MakeSpan(output.data<double>(), output.shape[0]);
auto size = std::min(outputSpan.size(), inputSpan.size());
REQUIRE( size > 0 );
@ -47,14 +47,14 @@ void testLowpass(const std::filesystem::path& inputNumpyFile, const std::filesys
OnePoleFilter filter { gain };
std::vector<Type> outputData (size);
filter.processLowpass(inputData, outputData);
filter.processLowpass(inputData, absl::MakeSpan(outputData));
REQUIRE( approxEqual(outputData, expectedData) );
filter.reset();
std::fill(outputData.begin(), outputData.end(), 0.0);
std::vector<Type> gains(size);
std::fill(gains.begin(), gains.end(), gain);
filter.processLowpassVariableGain(inputData, outputData, gains);
filter.processLowpassVariableGain(inputData, absl::MakeSpan(outputData), gains);
REQUIRE( approxEqual(outputData, expectedData) );
}
@ -63,11 +63,11 @@ void testHighpass(const std::filesystem::path& inputNumpyFile, const std::filesy
{
const auto input = cnpy::npy_load(inputNumpyFile.string());
REQUIRE( input.word_size == 8 );
const auto inputSpan = gsl::make_span(input.data<double>(), input.shape[0]);
const auto inputSpan = absl::MakeSpan(input.data<double>(), input.shape[0]);
const auto output = cnpy::npy_load(outputNumpyFile.string());
REQUIRE( output.word_size == 8 );
const auto outputSpan = gsl::make_span(output.data<double>(), output.shape[0]);
const auto outputSpan = absl::MakeSpan(output.data<double>(), output.shape[0]);
auto size = std::min(outputSpan.size(), inputSpan.size());
REQUIRE( size > 0 );
@ -82,14 +82,14 @@ void testHighpass(const std::filesystem::path& inputNumpyFile, const std::filesy
OnePoleFilter filter { gain };
std::vector<Type> outputData (size);
filter.processHighpass(inputData, outputData);
filter.processHighpass(inputData, absl::MakeSpan(outputData));
REQUIRE( approxEqual(outputData, expectedData) );
filter.reset();
std::fill(outputData.begin(), outputData.end(), 0.0);
std::vector<Type> gains(size);
std::fill(gains.begin(), gains.end(), gain);
filter.processHighpassVariableGain(inputData, outputData, gains);
filter.processHighpassVariableGain(inputData, absl::MakeSpan(outputData), gains);
REQUIRE( approxEqual(outputData, expectedData) );
}

View file

@ -1,5 +1,6 @@
#include "catch2/catch.hpp"
#include "../sources/SIMDHelpers.h"
#include <array>
#include <algorithm>
using namespace Catch::literals;
@ -12,7 +13,7 @@ TEST_CASE("[Helpers] fill() - Manual buffer")
{
std::vector<float> buffer (5);
std::vector<float> expected { fillValue, fillValue, fillValue, fillValue, fillValue };
fill<float, false>(buffer, fillValue);
fill<float, false>(absl::MakeSpan(buffer), fillValue);
REQUIRE(buffer == expected);
}
@ -22,7 +23,7 @@ TEST_CASE("[Helpers] fill() - Small buffer")
std::vector<float> expected (smallBufferSize);
std::fill(expected.begin(), expected.end(), fillValue);
fill<float, false>(buffer, fillValue);
fill<float, false>(absl::MakeSpan(buffer), fillValue);
REQUIRE(buffer == expected);
}
@ -32,7 +33,7 @@ TEST_CASE("[Helpers] fill() - Big buffer")
std::vector<float> expected (bigBufferSize);
std::fill(expected.begin(), expected.end(), fillValue);
fill<float, false>(buffer, fillValue);
fill<float, false>(absl::MakeSpan(buffer), fillValue);
REQUIRE(buffer == expected);
}
@ -42,7 +43,7 @@ TEST_CASE("[Helpers] fill() - Small buffer -- SIMD")
std::vector<float> expected (smallBufferSize);
std::fill(expected.begin(), expected.end(), fillValue);
fill<float, true>(buffer, fillValue);
fill<float, true>(absl::MakeSpan(buffer), fillValue);
REQUIRE(buffer == expected);
}
@ -52,7 +53,7 @@ TEST_CASE("[Helpers] fill() - Big buffer -- SIMD")
std::vector<float> expected (bigBufferSize);
std::fill(expected.begin(), expected.end(), fillValue);
fill<float, true>(buffer, fillValue);
fill<float, true>(absl::MakeSpan(buffer), fillValue);
REQUIRE(buffer == expected);
}
@ -62,7 +63,7 @@ TEST_CASE("[Helpers] fill() - Small buffer -- doubles")
std::vector<double> expected (smallBufferSize);
std::fill(expected.begin(), expected.end(), fillValue);
fill<double, false>(buffer, fillValue);
fill<double, false>(absl::MakeSpan(buffer), fillValue);
REQUIRE(buffer == expected);
}
@ -72,19 +73,19 @@ TEST_CASE("[Helpers] fill() - Big buffer -- doubles")
std::vector<double> expected (bigBufferSize);
std::fill(expected.begin(), expected.end(), fillValue);
fill<double, false>(buffer, fillValue);
fill<double, false>(absl::MakeSpan(buffer), fillValue);
REQUIRE(buffer == expected);
}
TEST_CASE("[Helpers] Interleaved read")
{
std::array<float, 16> input = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f};
std::array<float, 16> expected = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f };
std::array<float, 8> leftOutput { 0.0f };
std::array<float, 8> rightOutput { 0.0f };
readInterleaved<float, false>(input, leftOutput, rightOutput);
std::array<float, 16> real { 0.0f };
std::array<float, 16> input { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f};
std::array<float, 16> expected { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f };
std::array<float, 8> leftOutput;
std::array<float, 8> rightOutput;
readInterleaved<float, false>(input, absl::MakeSpan(leftOutput), absl::MakeSpan(rightOutput));
std::array<float, 16> real;
auto realIdx = 0;
for (auto value: leftOutput)
@ -96,12 +97,12 @@ TEST_CASE("[Helpers] Interleaved read")
TEST_CASE("[Helpers] Interleaved read unaligned end")
{
std::array<float, 20> input = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f, 8.0f, 18.0f, 9.0f, 19.0f};
std::array<float, 20> expected = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f, 18.0f, 19.0f};
std::array<float, 10> leftOutput { 0.0f };
std::array<float, 10> rightOutput { 0.0f };
readInterleaved<float, false>(input, leftOutput, rightOutput);
std::array<float, 20> real { 0.0f };
std::array<float, 20> input { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f, 8.0f, 18.0f, 9.0f, 19.0f};
std::array<float, 20> expected { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f, 18.0f, 19.0f};
std::array<float, 10> leftOutput;
std::array<float, 10> rightOutput;
readInterleaved<float, false>(input, absl::MakeSpan(leftOutput), absl::MakeSpan(rightOutput));
std::array<float, 20> real;
auto realIdx = 0;
for (auto value: leftOutput)
@ -113,12 +114,12 @@ TEST_CASE("[Helpers] Interleaved read unaligned end")
TEST_CASE("[Helpers] Small interleaved read unaligned end")
{
std::array<float, 6> input = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f};
std::array<float, 6> expected = { 0.0f, 1.0f, 2.0f, 10.0f, 11.0f, 12.0f};
std::array<float, 3> leftOutput { 0.0f };
std::array<float, 3> rightOutput { 0.0f };
readInterleaved<float, false>(input, leftOutput, rightOutput);
std::array<float, 6> real { 0.0f };
std::array<float, 6> input { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f};
std::array<float, 6> expected { 0.0f, 1.0f, 2.0f, 10.0f, 11.0f, 12.0f};
std::array<float, 3> leftOutput;
std::array<float, 3> rightOutput;
readInterleaved<float, false>(input, absl::MakeSpan(leftOutput), absl::MakeSpan(rightOutput));
std::array<float, 6> real;
auto realIdx = 0;
for (auto value: leftOutput)
@ -132,10 +133,10 @@ TEST_CASE("[Helpers] Interleaved read -- SIMD")
{
std::array<float, 16> input = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f};
std::array<float, 16> expected = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f };
std::array<float, 8> leftOutput { 0.0f };
std::array<float, 8> rightOutput { 0.0f };
readInterleaved<float, true>(input, leftOutput, rightOutput);
std::array<float, 16> real { 0.0f };
std::array<float, 8> leftOutput;
std::array<float, 8> rightOutput;
readInterleaved<float, true>(input, absl::MakeSpan(leftOutput), absl::MakeSpan(rightOutput));
std::array<float, 16> real;
auto realIdx = 0;
for (auto value: leftOutput)
@ -149,10 +150,10 @@ TEST_CASE("[Helpers] Interleaved read unaligned end -- SIMD")
{
std::array<float, 20> input = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f, 8.0f, 18.0f, 9.0f, 19.0f};
std::array<float, 20> expected = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f, 18.0f, 19.0f};
std::array<float, 10> leftOutput { 0.0f };
std::array<float, 10> rightOutput { 0.0f };
readInterleaved<float, true>(input, leftOutput, rightOutput);
std::array<float, 20> real { 0.0f };
std::array<float, 10> leftOutput;
std::array<float, 10> rightOutput;
readInterleaved<float, true>(input, absl::MakeSpan(leftOutput), absl::MakeSpan(rightOutput));
std::array<float, 20> real;
auto realIdx = 0;
for (auto value: leftOutput)
@ -164,12 +165,12 @@ TEST_CASE("[Helpers] Interleaved read unaligned end -- SIMD")
TEST_CASE("[Helpers] Small interleaved read unaligned end -- SIMD")
{
std::array<float, 6> input = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f};
std::array<float, 6> expected = { 0.0f, 1.0f, 2.0f, 10.0f, 11.0f, 12.0f};
std::array<float, 3> leftOutput { 0.0f };
std::array<float, 3> rightOutput { 0.0f };
readInterleaved<float, true>(input, leftOutput, rightOutput);
std::array<float, 6> real { 0.0f };
std::array<float, 6> input { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f};
std::array<float, 6> expected { 0.0f, 1.0f, 2.0f, 10.0f, 11.0f, 12.0f};
std::array<float, 3> leftOutput;
std::array<float, 3> rightOutput;
readInterleaved<float, true>(input, absl::MakeSpan(leftOutput), absl::MakeSpan(rightOutput));
std::array<float, 6> real;
auto realIdx = 0;
for (auto value: leftOutput)
@ -187,8 +188,8 @@ TEST_CASE("[Helpers] Interleaved read SIMD vs Scalar")
std::array<float, medBufferSize> leftOutputSIMD;
std::array<float, medBufferSize> rightOutputSIMD;
std::iota(input.begin(), input.end(), 0.0f);
readInterleaved<float, false>(input, leftOutputScalar, rightOutputScalar);
readInterleaved<float, true>(input, leftOutputSIMD, rightOutputSIMD);
readInterleaved<float, false>(input, absl::MakeSpan(leftOutputScalar), absl::MakeSpan(rightOutputScalar));
readInterleaved<float, true>(input, absl::MakeSpan(leftOutputSIMD), absl::MakeSpan(rightOutputSIMD));
REQUIRE( leftOutputScalar == leftOutputSIMD );
REQUIRE( rightOutputScalar == rightOutputSIMD );
}
@ -198,8 +199,8 @@ TEST_CASE("[Helpers] Interleaved write")
std::array<float, 8> leftInput { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, };
std::array<float, 8> rightInput { 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f };
std::array<float, 16> output;
std::array<float, 16> expected = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f};
writeInterleaved<float, false>(leftInput, rightInput, output);
std::array<float, 16> expected { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f};
writeInterleaved<float, false>(leftInput, rightInput, absl::MakeSpan(output));
REQUIRE( output == expected );
}
@ -208,8 +209,8 @@ TEST_CASE("[Helpers] Interleaved write unaligned end")
std::array<float, 10> leftInput { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f};
std::array<float, 10> rightInput { 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f, 18.0f, 19.0f };
std::array<float, 20> output;
std::array<float, 20> expected = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f, 8.0f, 18.0f, 9.0f, 19.0f};
writeInterleaved<float, false>(leftInput, rightInput, output);
std::array<float, 20> expected { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f, 8.0f, 18.0f, 9.0f, 19.0f};
writeInterleaved<float, false>(leftInput, rightInput, absl::MakeSpan(output));
REQUIRE( output == expected );
}
@ -218,8 +219,8 @@ TEST_CASE("[Helpers] Small interleaved write unaligned end")
std::array<float, 3> leftInput { 0.0f, 1.0f, 2.0f};
std::array<float, 3> rightInput { 10.0f, 11.0f, 12.0f };
std::array<float, 6> output;
std::array<float, 6> expected = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f};
writeInterleaved<float, false>(leftInput, rightInput, output);
std::array<float, 6> expected { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f};
writeInterleaved<float, false>(leftInput, rightInput, absl::MakeSpan(output));
REQUIRE( output == expected );
}
@ -228,8 +229,8 @@ TEST_CASE("[Helpers] Interleaved write -- SIMD")
std::array<float, 8> leftInput { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, };
std::array<float, 8> rightInput { 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f };
std::array<float, 16> output;
std::array<float, 16> expected = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f};
writeInterleaved<float, true>(leftInput, rightInput, output);
std::array<float, 16> expected { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f};
writeInterleaved<float, true>(leftInput, rightInput, absl::MakeSpan(output));
REQUIRE( output == expected );
}
@ -239,7 +240,7 @@ TEST_CASE("[Helpers] Interleaved write unaligned end -- SIMD")
std::array<float, 10> rightInput { 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f, 18.0f, 19.0f };
std::array<float, 20> output;
std::array<float, 20> expected = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f, 8.0f, 18.0f, 9.0f, 19.0f};
writeInterleaved<float, true>(leftInput, rightInput, output);
writeInterleaved<float, true>(leftInput, rightInput, absl::MakeSpan(output));
REQUIRE( output == expected );
}
@ -248,8 +249,8 @@ TEST_CASE("[Helpers] Small interleaved write unaligned end -- SIMD")
std::array<float, 3> leftInput { 0.0f, 1.0f, 2.0f};
std::array<float, 3> rightInput { 10.0f, 11.0f, 12.0f };
std::array<float, 6> output;
std::array<float, 6> expected = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f};
writeInterleaved<float, true>(leftInput, rightInput, output);
std::array<float, 6> expected { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f};
writeInterleaved<float, true>(leftInput, rightInput, absl::MakeSpan(output));
REQUIRE( output == expected );
}
@ -261,7 +262,7 @@ TEST_CASE("[Helpers] Interleaved write SIMD vs Scalar")
std::array<float, medBufferSize * 2> outputSIMD;
std::iota(leftInput.begin(), leftInput.end(), 0.0f);
std::iota(rightInput.begin(), rightInput.end(), medBufferSize);
writeInterleaved<float, false>(leftInput, rightInput, outputScalar);
writeInterleaved<float, true>(leftInput, rightInput, outputSIMD);
writeInterleaved<float, false>(leftInput, rightInput, absl::MakeSpan(outputScalar));
writeInterleaved<float, true>(leftInput, rightInput, absl::MakeSpan(outputSIMD));
REQUIRE( outputScalar == outputSIMD );
}

View file

@ -163,7 +163,7 @@ TEST_CASE("[StereoBuffer] Interleaved write -- Scalar")
std::array<float, 20> input = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f, 18.0f, 19.0f};
std::array<float, 20> output { 0.0f };
buffer.readInterleaved(input);
buffer.writeInterleaved(output);
buffer.writeInterleaved(absl::MakeSpan(output));
REQUIRE( output == input );
}
@ -173,7 +173,7 @@ TEST_CASE("[StereoBuffer] Interleaved write -- SIMD")
std::array<float, 20> input = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f, 3.0f, 13.0f, 4.0f, 14.0f, 5.0f, 15.0f, 6.0f, 16.0f, 7.0f, 17.0f, 8.0f, 18.0f, 9.0f, 19.0f};
std::array<float, 20> output { 0.0f };
buffer.readInterleaved(input);
buffer.writeInterleaved(output);
buffer.writeInterleaved(absl::MakeSpan(output));
REQUIRE( output == input );
}
@ -183,6 +183,6 @@ TEST_CASE("[StereoBuffer] Small interleaved write -- SIMD")
std::array<float, 6> input = { 0.0f, 10.0f, 1.0f, 11.0f, 2.0f, 12.0f};
std::array<float, 6> output { 0.0f };
buffer.readInterleaved(input);
buffer.writeInterleaved(output);
buffer.writeInterleaved(absl::MakeSpan(output));
REQUIRE( output == input );
}