view foosdk/sdk/pfc/audio_math.cpp @ 1:20d02a178406 default tip

*: check in everything else yay
author Paper <paper@tflc.us>
date Mon, 05 Jan 2026 02:15:46 -0500
parents
children
line wrap: on
line source

#include "pfc-lite.h"
#include "audio_sample.h"
#include "primitives.h"
#include "cpuid.h"


#if (defined(_M_IX86_FP) && _M_IX86_FP >= 2) || (defined(_M_X64) && !defined(_M_ARM64EC)) || defined(__x86_64__) || defined(__SSE2__)
#define AUDIO_MATH_SSE
#include <xmmintrin.h>
#include <tmmintrin.h> // _mm_shuffle_epi8
#include <smmintrin.h> // _mm_blend_epi16

#ifndef _mm_loadu_si32
#define _mm_loadu_si32(p) _mm_cvtsi32_si128(*(unsigned int const*)(p))
#endif
#ifndef _mm_storeu_si32
#define _mm_storeu_si32(p, a) (void)(*(int*)(p) = _mm_cvtsi128_si32((a)))
#endif

#ifdef __AVX__
#define allowAVX 1
#define haveAVX 1
#elif PFC_HAVE_CPUID
#define allowAVX 1
static const bool haveAVX = pfc::query_cpu_feature_set(pfc::CPU_HAVE_AVX);
#else
#define allowAVX 0
#define haveAVX 0
#endif

#ifdef __SSE4_1__
#define haveSSE41 true
#elif PFC_HAVE_CPUID
static const bool haveSSE41 = pfc::query_cpu_feature_set(pfc::CPU_HAVE_SSE41);
#else
#define haveSSE41 false
#endif

#if allowAVX
#include <immintrin.h> // _mm256_set1_pd
#endif

#endif // end SSE

#if defined( __aarch64__ ) || defined( _M_ARM64) || defined( _M_ARM64EC )
#define AUDIO_MATH_ARM64
#endif

#if defined( AUDIO_MATH_ARM64 ) || defined( __ARM_NEON__ )
#define AUDIO_MATH_NEON
#include <arm_neon.h>

// No vcvtnq_s32_f32 on ARM32, use vcvtq_s32_f32, close enough
#ifdef AUDIO_MATH_ARM64
#define vcvtnq_s32_f32_wrap vcvtnq_s32_f32
#else
#define vcvtnq_s32_f32_wrap vcvtq_s32_f32
#endif

#endif


#if defined( AUDIO_MATH_ARM64 ) && !defined( __ANDROID__ )
// Don't do Neon float64 on Android, crashes clang from NDK 25
#define AUDIO_MATH_NEON_FLOAT64
#endif

template<typename float_t> inline static float_t noopt_calculate_peak(const float_t *p_src, t_size p_num)
{
	float_t peak = 0;
	t_size num = p_num;
	for(;num;num--)
	{
		float_t temp = (float_t)fabs(*(p_src++));
        peak = fmax(peak, temp);
	}
	return peak;
}


template<typename float_t>
inline static void noopt_convert_to_32bit(const float_t* p_source,t_size p_count,t_int32 * p_output, float_t p_scale)
{
	t_size num = p_count;
	for(;num;--num)
	{
		t_int64 val = pfc::audio_math::rint64( *(p_source++) * p_scale );
		if (val < INT32_MIN) val = INT32_MIN;
		else if (val > INT32_MAX) val = INT32_MAX;
		*(p_output++) = (t_int32) val;
	}
}

template<typename float_t>
inline static void noopt_convert_to_16bit(const float_t* p_source,t_size p_count,t_int16 * p_output, float_t p_scale) {
	for(t_size n=0;n<p_count;n++) {
		*(p_output++) = (t_int16) pfc::clip_t<int32_t>(pfc::audio_math::rint32(*(p_source++)*p_scale),INT16_MIN,INT16_MAX);
	}
}

template<typename float_t>
inline static void noopt_convert_from_int16(const t_int16 * __restrict p_source,t_size p_count, float_t* __restrict p_output, float_t p_scale)
{
	t_size num = p_count;
	for(;num;num--)
		*(p_output++) = (float_t)*(p_source++) * p_scale;
}



template<typename float_t>
inline static void noopt_convert_from_int32(const t_int32 * __restrict p_source,t_size p_count, float_t* __restrict p_output, float_t p_scale)
{
	t_size num = p_count;
	for(;num;num--)
		*(p_output++) = (float_t)( * (p_source++) * p_scale );
}

template<typename in_t, typename out_t, typename scale_t>
inline static void noopt_scale(const in_t * p_source,size_t p_count,out_t * p_output,scale_t p_scale)
{
	for(t_size n=0;n<p_count;n++)
		p_output[n] = (out_t)(p_source[n] * p_scale);
}
template<typename in_t, typename out_t>
inline static void noopt_convert(const in_t* in, out_t* out, size_t count) {
    for (size_t walk = 0; walk < count; ++walk) out[walk] = (out_t)in[walk];
}

#ifdef AUDIO_MATH_NEON

#ifdef AUDIO_MATH_ARM64
#define _vmaxvq_f32_wrap vmaxvq_f32
#else
inline float _vmaxvq_f32_wrap( float32x4_t arg ) {
    return pfc::max_t<float>( pfc::max_t<float>(arg[0], arg[1]), pfc::max_t<float>(arg[2], arg[3]) );
}
#endif

inline static float neon_calculate_peak( const float * p_source, size_t p_count ) {
    size_t num = p_count / 8;
    float32x4_t ret1 = {}, ret2 = {};
    for(;num;--num) {
        float32x4_t f32lo = vld1q_f32( p_source );
        float32x4_t f32hi = vld1q_f32( p_source + 4 );
        p_source += 8;
        ret1 = vmaxq_f32(ret1, vabsq_f32(f32lo));
        ret2 = vmaxq_f32(ret2, vabsq_f32(f32hi));
    }
    
    float ret = _vmaxvq_f32_wrap(vmaxq_f32( ret1, ret2 ));
    
    size_t rem = p_count % 8;
    if ( rem != 0 ) {
        float v = noopt_calculate_peak( p_source, p_count % 8);
        if (v > ret) ret = v;
    }

    return ret;
}

inline static void neon_scale(const float * p_source,size_t p_count, float * p_output,float p_scale) {
    size_t num = p_count / 8;
    for(;num;--num) {
        float32x4_t lo = vld1q_f32( p_source );
        float32x4_t hi = vld1q_f32( p_source + 4 );

        lo = vmulq_n_f32( lo, p_scale);
        hi = vmulq_n_f32( hi, p_scale);
        
        vst1q_f32( p_output, lo );
        vst1q_f32( p_output+4, hi );
        
        p_source += 8;
        p_output += 8;
    }
    
    noopt_scale( p_source, p_count % 8, p_output, p_scale);
}
inline static void neon_convert_to_int32(const float * __restrict p_source,t_size p_count, int32_t * __restrict p_output,float p_scale)
{
    size_t num = p_count / 8;
    for(;num;--num) {
        float32x4_t f32lo = vld1q_f32( p_source );
        float32x4_t f32hi = vld1q_f32( p_source + 4 );
        
        int32x4_t lo = vcvtnq_s32_f32_wrap( vmulq_n_f32(f32lo, p_scale) );
        int32x4_t hi = vcvtnq_s32_f32_wrap( vmulq_n_f32(f32hi, p_scale) );
        
        vst1q_s32(p_output, lo);
        vst1q_s32(p_output+4, hi);
        
        p_source += 8;
        p_output += 8;
        
    }
    
    noopt_convert_to_32bit(p_source, p_count % 8, p_output, p_scale);
}

inline static void neon_convert_from_int32(const int32_t * __restrict p_source,t_size p_count, float * __restrict p_output,float p_scale)
{
    size_t num = p_count / 8;
    size_t rem = p_count % 8;
    for(;num;num--) {
        int32x4_t i32lo = vld1q_s32( p_source );
        int32x4_t i32hi = vld1q_s32( p_source + 4 );
        float32x4_t f32vl = vcvtq_f32_s32(i32lo);
        float32x4_t f32vh = vcvtq_f32_s32(i32hi);

        vst1q_f32(&p_output[0], vmulq_n_f32(f32vl, p_scale));
        vst1q_f32(&p_output[4], vmulq_n_f32(f32vh, p_scale));
        
        p_source += 8;
        p_output += 8;
        
    }
    
    noopt_convert_from_int32( p_source, rem, p_output, p_scale );
}

inline static void neon_convert_to_int16(const float* __restrict p_source,t_size p_count, int16_t * __restrict p_output,float p_scale)
{
    size_t num = p_count / 8;
    size_t rem = p_count % 8;
    for(;num;--num) {
        float32x4_t f32lo = vld1q_f32( p_source );
        float32x4_t f32hi = vld1q_f32( p_source + 4);
        
        int32x4_t lo = vcvtnq_s32_f32_wrap( vmulq_n_f32(f32lo, p_scale) );
        int32x4_t hi = vcvtnq_s32_f32_wrap( vmulq_n_f32(f32hi, p_scale) );
        
        vst1q_s16(&p_output[0], vcombine_s16( vqmovn_s32( lo ), vqmovn_s32( hi ) ) );
        
        p_source += 8;
        p_output += 8;
        
    }
    
    noopt_convert_to_16bit(p_source, rem, p_output, p_scale);
    
}
inline static void neon_convert_from_int16(const t_int16 * __restrict p_source,t_size p_count, float * __restrict p_output,float p_scale)
{
    size_t num = p_count / 8;
    size_t rem = p_count % 8;
    for(;num;num--) {
        auto i16lo = vld1_s16(p_source);
        auto i16hi = vld1_s16(p_source + 4);

        float32x4_t f32vl = vcvtq_f32_s32(vmovl_s16 (i16lo));
        float32x4_t f32vh = vcvtq_f32_s32(vmovl_s16 (i16hi));

        vst1q_f32(&p_output[0], vmulq_n_f32(f32vl, p_scale));
        vst1q_f32(&p_output[4], vmulq_n_f32(f32vh, p_scale));

        p_source += 8;
        p_output += 8;

    }
    
    noopt_convert_from_int16( p_source, rem, p_output, p_scale );
}
#ifdef AUDIO_MATH_NEON_FLOAT64
inline static void neon_convert_to_int16(const double* __restrict p_source, t_size p_count, int16_t* __restrict p_output, double p_scale)
{
    size_t num = p_count / 4;
    size_t rem = p_count % 4;
    for (; num; --num) {
        float64x2_t f64lo = vld1q_f64(p_source);
        float64x2_t f64hi = vld1q_f64(p_source + 2);

        f64lo = vmulq_n_f64(f64lo, p_scale);
        f64hi = vmulq_n_f64(f64hi, p_scale);

        int64x2_t lo64 = vcvtnq_s64_f64(f64lo);
        int64x2_t hi64 = vcvtnq_s64_f64(f64hi);
        
        int32x4_t v32 = vcombine_s32(vqmovn_s64(lo64), vqmovn_s64(hi64));


        vst1_s16(&p_output[0], vqmovn_s32(v32));

        p_source += 4;
        p_output += 4;

    }

    noopt_convert_to_16bit(p_source, rem, p_output, p_scale);
}

inline static void neon_convert_from_int16(const t_int16* __restrict p_source, t_size p_count, double* __restrict p_output, double p_scale)
{
    size_t num = p_count / 4;
    size_t rem = p_count % 4;
    for (; num; num--) {
        int32x4_t i32 = vmovl_s16(vld1_s16(p_source));
        
        int64x2_t lo64 = vmovl_s32( vget_low_s32(i32) );
        int64x2_t hi64 = vmovl_s32(vget_high_s32(i32));

        float64x2_t f64vl = vcvtq_f64_s64(lo64);
        float64x2_t f64vh = vcvtq_f64_s64(hi64);

        vst1q_f64(&p_output[0], vmulq_n_f64(f64vl, p_scale));
        vst1q_f64(&p_output[2], vmulq_n_f64(f64vh, p_scale));

        p_source += 4;
        p_output += 4;

    }

    noopt_convert_from_int16(p_source, rem, p_output, p_scale);
}
#endif // AUDIO_MATH_NEON_FLOAT64

#endif // AUDIO_MATH_NEON

#if defined(AUDIO_MATH_SSE)

inline void convert_to_32bit_sse2(const float* p_src, size_t numTotal, t_int32* p_dst, float p_mul)
{

    // Implementation notes
    // There doesn't seem to be a nice and tidy way to convert float to int32 with graceful clipping to INT32_MIN .. INT32_MAX range.
    // While low clipping at INT32_MIN can be accomplished with _mm_max_ps(), high clipping needs float compare THEN substitute bad int with INT32_MAX.
    // The best we could do with _mm_min_ps() would result with high clipping at 0x7FFFFF80 instead of 0x7FFFFFFF (INT32_MAX).
    // We store masks from float compare and fix ints according to the mask later.

    __m128 mul = _mm_set1_ps(p_mul);
    __m128 loF = _mm_set1_ps((float)INT32_MIN);
    __m128 hiF = _mm_set1_ps((float)INT32_MAX);
    // __m128i loI = _mm_set1_epi32(INT32_MIN);
    __m128i hiI = _mm_set1_epi32(INT32_MAX);

    size_t num = numTotal / 4;
    size_t rem = numTotal % 4;
    for (; num; --num) {
        __m128 s = _mm_mul_ps(mul, _mm_loadu_ps(p_src)); p_src += 4;

        s = _mm_max_ps(s, loF);

        // __m128i maskLo = _mm_castps_si128(_mm_cmple_ps(s, loF));
        __m128i maskHi = _mm_castps_si128(_mm_cmpge_ps(s, hiF));

        __m128i i = _mm_cvtps_epi32(s);

        // i = _mm_or_si128(_mm_andnot_si128(maskLo, i), _mm_and_si128(loI, maskLo));
        i = _mm_or_si128(_mm_andnot_si128(maskHi, i), _mm_and_si128(hiI, maskHi));

        _mm_storeu_si128((__m128i*) p_dst, i); p_dst += 4;
    }

    for (; rem; --rem) {
        __m128 s = _mm_mul_ss(_mm_load_ss(p_src++), mul);
        s = _mm_max_ss(s, loF);

        // __m128i maskLo = _mm_castps_si128( _mm_cmple_ss(s, loF) );
        __m128i maskHi = _mm_castps_si128(_mm_cmpge_ss(s, hiF));

        __m128i i = _mm_cvtps_epi32(s); // not ss

        // i = _mm_or_si128(_mm_andnot_si128(maskLo, i), _mm_and_si128(loI, maskLo));
        i = _mm_or_si128(_mm_andnot_si128(maskHi, i), _mm_and_si128(hiI, maskHi));

        _mm_storeu_si32(p_dst++, i);
    }
}

inline void convert_to_32bit_sse2(const double* p_src, size_t numTotal, t_int32* p_dst, double p_mul)
{
    auto mul = _mm_set1_pd(p_mul);
    auto loF = _mm_set1_pd(INT32_MIN);
    auto hiF = _mm_set1_pd(INT32_MAX);

    size_t num = numTotal / 4;
    size_t rem = numTotal % 4;
    for (; num; --num) {
        auto v1 = _mm_loadu_pd(p_src);
        auto v2 = _mm_loadu_pd(p_src + 2);
        p_src += 4;

        v1 = _mm_mul_pd(v1, mul); v2 = _mm_mul_pd(v2, mul);

        v1 = _mm_max_pd(v1, loF); v2 = _mm_max_pd(v2, loF);
        v1 = _mm_min_pd(v1, hiF); v2 = _mm_min_pd(v2, hiF);

        auto i1 = _mm_cvtpd_epi32(v1), i2 = _mm_cvtpd_epi32(v2);
        

        _mm_storeu_si128((__m128i*) p_dst, _mm_unpacklo_epi64(i1, i2)); p_dst += 4;
    }

    for (; rem; --rem) {
        auto s = _mm_mul_sd(_mm_load_sd(p_src++), mul);
        s = _mm_max_sd(s, loF); s = _mm_min_sd(s, hiF);
        * p_dst++ = _mm_cvtsd_si32(s);
    }
}

inline void convert_from_int16_sse2(const t_int16 * p_source,t_size p_count,float * p_output,float p_scale)
{
    while(!pfc::is_ptr_aligned_t<16>(p_output) && p_count) {
        *(p_output++) = (float)*(p_source++) * p_scale;
        p_count--;
    }

    {
        __m128 mul = _mm_set1_ps(p_scale);
        __m128i nulls = _mm_setzero_si128();
        __m128i delta1 = _mm_set1_epi16((int16_t)0x8000);
        __m128i delta2 = _mm_set1_epi32((int32_t)0x8000);

        for(t_size loop = p_count >> 3;loop;--loop) {
            __m128i source, temp1, temp2; __m128 float1, float2;
            source = _mm_loadu_si128((__m128i*)p_source);
            source = _mm_xor_si128(source,delta1);
            temp1 = _mm_unpacklo_epi16(source,nulls);
            temp2 = _mm_unpackhi_epi16(source,nulls);
            temp1 = _mm_sub_epi32(temp1,delta2);
            temp2 = _mm_sub_epi32(temp2,delta2);
            p_source += 8;
            float1 = _mm_cvtepi32_ps(temp1);
            float2 = _mm_cvtepi32_ps(temp2);
            float1 = _mm_mul_ps(float1,mul);
            float2 = _mm_mul_ps(float2,mul);
            _mm_store_ps(p_output,float1);
            _mm_store_ps(p_output+4,float2);
            p_output += 8;
        }
        
        p_count &= 7;
    }

    while(p_count) {
        *(p_output++) = (float)*(p_source++) * p_scale;
        p_count--;
    }
}

inline static void convert_to_16bit_sse2(const float * p_source,t_size p_count,t_int16 * p_output,float p_scale)
{
    size_t num = p_count/8;
    size_t rem = p_count%8;
    __m128 mul = _mm_set1_ps(p_scale);
    for(;num;--num)
    {
        __m128 temp1,temp2; __m128i itemp1, itemp2;
        temp1 = _mm_loadu_ps(p_source);
        temp2 = _mm_loadu_ps(p_source+4);
        temp1 = _mm_mul_ps(temp1,mul);
        temp2 = _mm_mul_ps(temp2,mul);
        p_source += 8;
        itemp1 = _mm_cvtps_epi32(temp1);
        itemp2 = _mm_cvtps_epi32(temp2);
        _mm_storeu_si128( (__m128i*)p_output, _mm_packs_epi32(itemp1, itemp2) );
        p_output += 8;
    }
    
    noopt_convert_to_16bit(p_source, rem, p_output, p_scale);
}

inline static void convert_to_16bit_sse2(const double* p_source, t_size p_count, t_int16* p_output, double p_scale)
{
    size_t num = p_count / 8;
    size_t rem = p_count % 8;
    __m128d mul = _mm_set1_pd(p_scale);
    for (; num; --num)
    {
        __m128d temp1, temp2, temp3, temp4; __m128i itemp1, itemp2;
        temp1 = _mm_loadu_pd(p_source);
        temp2 = _mm_loadu_pd(p_source + 2);
        temp3 = _mm_loadu_pd(p_source + 4);
        temp4 = _mm_loadu_pd(p_source + 6);

        p_source += 8;

        temp1 = _mm_mul_pd(temp1, mul);
        temp2 = _mm_mul_pd(temp2, mul);
        temp3 = _mm_mul_pd(temp3, mul);
        temp4 = _mm_mul_pd(temp4, mul);

        
        itemp1 = _mm_unpacklo_epi64(_mm_cvtpd_epi32(temp1), _mm_cvtpd_epi32(temp2));
        itemp2 = _mm_unpacklo_epi64(_mm_cvtpd_epi32(temp3), _mm_cvtpd_epi32(temp4));

        _mm_storeu_si128((__m128i*)p_output, _mm_packs_epi32(itemp1, itemp2));
        p_output += 8;
    }

    noopt_convert_to_16bit(p_source, rem, p_output, p_scale);
}
#if allowAVX
inline static void avx_convert_to_16bit(const double* p_source, size_t p_count, int16_t* p_output, double p_scale) {
    size_t num = p_count / 8;
    size_t rem = p_count % 8;
    auto mul = _mm256_set1_pd(p_scale);
    for (; num; --num)
    {
        auto temp1 = _mm256_loadu_pd(p_source);
        auto temp2 = _mm256_loadu_pd(p_source + 4);

        p_source += 8;

        temp1 = _mm256_mul_pd(temp1, mul);
        temp2 = _mm256_mul_pd(temp2, mul);

        auto itemp1 = _mm256_cvtpd_epi32(temp1);
        auto itemp2 = _mm256_cvtpd_epi32(temp2);

        _mm_storeu_si128((__m128i*)p_output, _mm_packs_epi32(itemp1, itemp2));
        p_output += 8;
    }

    noopt_convert_to_16bit(p_source, rem, p_output, p_scale);
}
#endif

inline float sse_calculate_peak( const float * src, size_t count ) {
    size_t num = count/8;
    size_t rem = count%8;
        
    __m128 mask = _mm_castsi128_ps(_mm_set1_epi32(0x7FFFFFFF));
    __m128 acc1 = _mm_setzero_ps(), acc2 = _mm_setzero_ps();
    
    for(;num;--num) {
        __m128 v1 = _mm_loadu_ps( src );
        __m128 v2 = _mm_loadu_ps( src + 4 );
        v1 = _mm_and_ps( v1, mask );
        v2 = _mm_and_ps( v2, mask );
        // Two acc channels so one _mm_max_ps doesn't block the other
        acc1 = _mm_max_ps( acc1, v1 );
        acc2 = _mm_max_ps( acc2, v2 );
        src += 8;
    }
    
    float ret;
    {
        float blah[4];
        _mm_storeu_ps(blah, _mm_max_ps( acc1, acc2 ));
        __m128 acc = _mm_load_ss( &blah[0] );
        acc = _mm_max_ss( acc, _mm_load_ss( &blah[1] ) );
        acc = _mm_max_ss( acc, _mm_load_ss( &blah[2] ) );
        acc = _mm_max_ss( acc, _mm_load_ss( &blah[3] ) );
        ret = _mm_cvtss_f32(acc);
    }
    if ( rem > 0 ) {
        __m128 acc = _mm_set_ss( ret );
        for( ;rem; --rem) {
            __m128 v = _mm_load_ss( src++ );
            v = _mm_and_ps( v, mask );
            acc = _mm_max_ss( acc, v );
        }
        ret = _mm_cvtss_f32(acc);
    }
    return ret;
}

#if allowAVX
inline double avx_calculate_peak(const double* src, size_t count) {
    size_t num = count / 8;
    size_t rem = count % 8;

    auto mask = _mm256_castsi256_pd(_mm256_set1_epi64x(0x7FFFFFFFFFFFFFFF));
    auto acc1 = _mm256_setzero_pd(), acc2 = _mm256_setzero_pd();

    for (; num; --num) {
        auto v1 = _mm256_loadu_pd(src);
        auto v2 = _mm256_loadu_pd(src + 4);

        v1 = _mm256_and_pd(v1, mask);
        v2 = _mm256_and_pd(v2, mask);

        acc1 = _mm256_max_pd(acc1, v1);
        acc2 = _mm256_max_pd(acc2, v2);

        src += 8;
    }

    __m128d acc;
    {
        acc1 = _mm256_max_pd(acc1, acc2);

        acc = _mm_max_pd(_mm256_extractf128_pd(acc1, 0), _mm256_extractf128_pd(acc1, 1));

        acc = _mm_max_sd(acc, _mm_shuffle_pd(acc, acc, _MM_SHUFFLE2(0, 1)));
    }

    if (rem > 0) {
        __m128d mask128 = _mm_castsi128_pd(_mm_set1_epi64x(0x7FFFFFFFFFFFFFFF));
        for (; rem; --rem) {
            __m128d v = _mm_load_sd(src++);
            v = _mm_and_pd(v, mask128);
            acc = _mm_max_sd(acc, v);
        }
    }
    return _mm_cvtsd_f64(acc);
}
#endif // allowAVX

inline double sse_calculate_peak(const double* src, size_t count) {
    size_t num = count / 4;
    size_t rem = count % 4;

    __m128d mask = _mm_castsi128_pd(_mm_set1_epi64x(0x7FFFFFFFFFFFFFFF));
    __m128d acc1 = _mm_setzero_pd(), acc2 = _mm_setzero_pd();

    for (; num; --num) {
        __m128d v1 = _mm_loadu_pd(src);
        __m128d v2 = _mm_loadu_pd(src + 2);
        v1 = _mm_and_pd(v1, mask);
        v2 = _mm_and_pd(v2, mask);
        // Two acc channels so one _mm_max_pd doesn't block the other
        acc1 = _mm_max_pd(acc1, v1);
        acc2 = _mm_max_pd(acc2, v2);
        src += 4;
    }

    {
        acc1 = _mm_max_pd(acc1, acc2);
        acc1 = _mm_max_sd(acc1, _mm_shuffle_pd(acc1, acc1, _MM_SHUFFLE2(0, 1)));
    }
    if (rem > 0) {
        for (; rem; --rem) {
            __m128d v = _mm_load_sd(src++);
            v = _mm_and_pd(v, mask);
            acc1 = _mm_max_sd(acc1, v);
        }
    }
    return _mm_cvtsd_f64(acc1);
}

inline void sse_convert_from_int32(const int32_t* source, size_t count, float* output, float scale) {
    __m128 mul = _mm_set1_ps(scale);
    for (size_t num = count/8; num; --num)
    {
        __m128i itemp1, itemp2; __m128 temp1, temp2;
        itemp1 = _mm_loadu_si128((__m128i*)source);
        itemp2 = _mm_loadu_si128((__m128i*)source + 1);
        temp1 = _mm_cvtepi32_ps(itemp1);
        temp2 = _mm_cvtepi32_ps(itemp2);
        source += 8;
        temp1 = _mm_mul_ps(temp1, mul);
        temp2 = _mm_mul_ps(temp2, mul);
        _mm_storeu_ps(output, temp1);
        _mm_storeu_ps(output + 4, temp2);
        output += 8;
    }
    for (size_t rem = count % 8; rem; --rem) {
        __m128i i = _mm_loadu_si32(source++);
        __m128 f = _mm_cvtepi32_ps(i);
        f = _mm_mul_ss(f, mul);
        _mm_store_ss(output++, f);
    }
}

inline void sse_convert_from_int32(const int32_t* source, size_t count, double* output, double scale) {
    auto mul = _mm_set1_pd(scale);
    for (size_t num = count / 8; num; --num)
    {
        auto itemp1 = _mm_loadu_si128((__m128i*)source);
        auto itemp2 = _mm_loadu_si128((__m128i*)source + 1);
        auto temp1 = _mm_cvtepi32_pd(itemp1);
        auto temp2 = _mm_cvtepi32_pd(_mm_shuffle_epi32(itemp1, _MM_SHUFFLE(1, 0, 3, 2)));
        auto temp3 = _mm_cvtepi32_pd(itemp2);
        auto temp4 = _mm_cvtepi32_pd(_mm_shuffle_epi32(itemp2, _MM_SHUFFLE(1, 0, 3, 2)));
        source += 8;
        temp1 = _mm_mul_pd(temp1, mul);
        temp2 = _mm_mul_pd(temp2, mul);
        temp3 = _mm_mul_pd(temp3, mul);
        temp4 = _mm_mul_pd(temp4, mul);
        _mm_storeu_pd(output, temp1);
        _mm_storeu_pd(output + 2, temp2);
        _mm_storeu_pd(output + 4, temp3);
        _mm_storeu_pd(output + 6, temp4);
        output += 8;
    }
    for (size_t rem = count % 8; rem; --rem) {
        __m128i i = _mm_loadu_si32(source++);
        auto f = _mm_cvtepi32_pd(i);
        f = _mm_mul_sd(f, mul);
        _mm_store_sd(output++, f);
    }
}
#if allowAVX
inline void convert_from_int16_avx(const t_int16* p_source, t_size p_count, double* p_output, double p_scale) {
    while (!pfc::is_ptr_aligned_t<32>(p_output) && p_count) {
        *(p_output++) = (double)*(p_source++) * p_scale;
        p_count--;
    }

    {
        __m256d muld = _mm256_set1_pd(p_scale);

        for (t_size loop = p_count >> 3; loop; --loop) {
            auto source = _mm_loadu_si128((__m128i*)p_source);
            auto temp1 = _mm_cvtepi16_epi32(source);
            auto temp2 = _mm_cvtepi16_epi32(_mm_shuffle_epi32(source, _MM_SHUFFLE(0, 0, 3, 2)));
            p_source += 8;

            auto double1 = _mm256_cvtepi32_pd(temp1);
            auto double2 = _mm256_cvtepi32_pd(temp2);

            double1 = _mm256_mul_pd(double1, muld);
            double2 = _mm256_mul_pd(double2, muld);

            _mm256_store_pd(p_output, double1);
            _mm256_store_pd(p_output+4, double2);

            p_output += 8;
        }

        p_count &= 7;
    }

    while (p_count) {
        *(p_output++) = (double)*(p_source++) * p_scale;
        p_count--;
    }

}
#endif // allowAVX

inline void convert_from_int16_sse2(const t_int16* p_source, t_size p_count, double * p_output, double p_scale)
{
    while (!pfc::is_ptr_aligned_t<16>(p_output) && p_count) {
        *(p_output++) = (double) * (p_source++) * p_scale;
        p_count--;
    }

    {
        __m128d muld = _mm_set1_pd(p_scale);
        __m128i nulls = _mm_setzero_si128();
        __m128i delta1 = _mm_set1_epi16((int16_t)0x8000);
        __m128i delta2 = _mm_set1_epi32((int32_t)0x8000);

        for (t_size loop = p_count >> 3; loop; --loop) {
            __m128i source, temp1, temp2; __m128d double1, double2, double3, double4;
            source = _mm_loadu_si128((__m128i*)p_source);
            source = _mm_xor_si128(source, delta1);
            temp1 = _mm_unpacklo_epi16(source, nulls);
            temp2 = _mm_unpackhi_epi16(source, nulls);
            temp1 = _mm_sub_epi32(temp1, delta2);
            temp2 = _mm_sub_epi32(temp2, delta2);
            p_source += 8;

            double1 = _mm_cvtepi32_pd(temp1);
            double2 = _mm_cvtepi32_pd(_mm_shuffle_epi32(temp1, _MM_SHUFFLE(3, 2, 3, 2)));
            double3 = _mm_cvtepi32_pd(temp2);
            double4 = _mm_cvtepi32_pd(_mm_shuffle_epi32(temp2, _MM_SHUFFLE(3, 2, 3, 2)));

            double1 = _mm_mul_pd(double1, muld);
            double2 = _mm_mul_pd(double2, muld);
            double3 = _mm_mul_pd(double3, muld);
            double4 = _mm_mul_pd(double4, muld);
            _mm_store_pd(p_output, double1);
            _mm_store_pd(p_output + 2, double2);
            _mm_store_pd(p_output + 4, double3);
            _mm_store_pd(p_output + 6, double4);

            p_output += 8;
        }

        p_count &= 7;
    }

    while (p_count) {
        *(p_output++) = (double) * (p_source++) * p_scale;
        p_count--;
    }
}

#endif

namespace pfc {
    void audio_math::scale(const float* p_source, size_t p_count, float* p_output, float p_scale) {
#if defined( AUDIO_MATH_NEON )
        neon_scale(p_source, p_count, p_output, p_scale);
#else
        noopt_scale(p_source, p_count, p_output, p_scale);
#endif
    }
    void audio_math::scale(const double* p_source, size_t p_count, double* p_output, double p_scale) {
        noopt_scale(p_source, p_count, p_output, p_scale);
    }

    void audio_math::convert_to_int16(const float* p_source, t_size p_count, t_int16* p_output, float p_scale)
    {
        float scale = (float)(p_scale * 0x8000);
#if defined(AUDIO_MATH_SSE)
        convert_to_16bit_sse2(p_source, p_count, p_output, scale);
#elif defined( AUDIO_MATH_NEON )
        neon_convert_to_int16(p_source, p_count, p_output, scale);
#else
        noopt_convert_to_16bit(p_source, p_count, p_output, scale);
#endif
    }
    void audio_math::convert_to_int16(const double* p_source, t_size p_count, t_int16* p_output, double p_scale)
    {
        double scale = (double)(p_scale * 0x8000);
#if defined(AUDIO_MATH_SSE)
#if allowAVX
        if (haveAVX) {
            avx_convert_to_16bit(p_source, p_count, p_output, scale);
        } else
#endif // allowAVX
        {
            convert_to_16bit_sse2(p_source, p_count, p_output, scale);
        }
#elif defined( AUDIO_MATH_NEON_FLOAT64 )
        neon_convert_to_int16(p_source, p_count, p_output, scale);
#else
        noopt_convert_to_16bit(p_source, p_count, p_output, scale);
#endif
    }

    void audio_math::convert_from_int16(const t_int16* p_source, t_size p_count, float* p_output, float p_scale)
    {
        float scale = (float)(p_scale / (double)0x8000);
#if defined(AUDIO_MATH_SSE)
        convert_from_int16_sse2(p_source, p_count, p_output, scale);
#elif defined( AUDIO_MATH_NEON )
        neon_convert_from_int16(p_source, p_count, p_output, scale);
#else
        noopt_convert_from_int16(p_source, p_count, p_output, scale);
#endif
    }

    void audio_math::convert_from_int16(const t_int16* p_source, t_size p_count, double* p_output, double p_scale)
    {
        double scale = (double)(p_scale / (double)0x8000);
#if defined(AUDIO_MATH_SSE)
#if allowAVX
        if (haveAVX) {
            convert_from_int16_avx(p_source, p_count, p_output, scale);
        } else
#endif
        {
            convert_from_int16_sse2(p_source, p_count, p_output, scale);
        }
#elif defined( AUDIO_MATH_NEON_FLOAT64 )
        neon_convert_from_int16(p_source, p_count, p_output, scale);
#else
        noopt_convert_from_int16(p_source, p_count, p_output, scale);
#endif
    }

    void audio_math::convert_to_int32(const float* p_source, t_size p_count, t_int32* p_output, float p_scale)
    {
        float scale = (float)(p_scale * 0x80000000ul);
#if defined(AUDIO_MATH_NEON)
        neon_convert_to_int32(p_source, p_count, p_output, scale);
#elif defined(AUDIO_MATH_SSE)
        convert_to_32bit_sse2(p_source, p_count, p_output, scale);
#else
        noopt_convert_to_32bit(p_source, p_count, p_output, scale);
#endif
    }

    void audio_math::convert_to_int32(const double* p_source, t_size p_count, t_int32* p_output, double p_scale)
    {
        double scale = (double)(p_scale * 0x80000000ul);
#if defined(AUDIO_MATH_SSE)
        convert_to_32bit_sse2(p_source, p_count, p_output, scale);
#else
        noopt_convert_to_32bit(p_source, p_count, p_output, scale);
#endif
    }

    void audio_math::convert_from_int32(const t_int32* p_source, t_size p_count, float* p_output, float p_scale)
    {
        float scale = (float)(p_scale / (double)0x80000000ul);
        // Note: speed difference here is marginal over compiler output as of Xcode 12
#if defined(AUDIO_MATH_NEON)
        neon_convert_from_int32(p_source, p_count, p_output, scale);
#elif defined(AUDIO_MATH_SSE)
        sse_convert_from_int32(p_source, p_count, p_output, scale);
#else
        noopt_convert_from_int32(p_source, p_count, p_output, scale);
#endif
    }

    void audio_math::convert_from_int32(const t_int32* p_source, t_size p_count, double* p_output, double p_scale)
    {
        double scale = (double)(p_scale / (double)0x80000000ul);
#if defined(AUDIO_MATH_SSE)
        sse_convert_from_int32(p_source, p_count, p_output, scale);
#else
        noopt_convert_from_int32(p_source, p_count, p_output, scale);
#endif
    }

    float audio_math::calculate_peak(const float * p_source, t_size p_count) {
#if defined(AUDIO_MATH_SSE)
        return sse_calculate_peak(p_source, p_count);
#elif defined(AUDIO_MATH_NEON)
        return neon_calculate_peak(p_source, p_count);
#else
        return noopt_calculate_peak(p_source, p_count);
#endif
    }
    double audio_math::calculate_peak(const double * p_source, t_size p_count) {
#if defined(AUDIO_MATH_SSE)
        // Note that avx_calculate_peak failed to score better than sse_calculate_peak
        return sse_calculate_peak(p_source, p_count);
#else
        return noopt_calculate_peak(p_source, p_count);
#endif
    }

    void audio_math::remove_denormals(float* p_buffer, t_size p_count) {
        t_uint32* ptr = reinterpret_cast<t_uint32*>(p_buffer);
        for (; p_count; p_count--)
        {
            t_uint32 t = *ptr;
            if ((t & 0x007FFFFF) && !(t & 0x7F800000)) *ptr = 0;
            ptr++;
        }
    }
    void audio_math::remove_denormals(double* p_buffer, t_size p_count) {
        t_uint64* ptr = reinterpret_cast<t_uint64*>(p_buffer);
        for (; p_count; p_count--)
        {
            t_uint64 t = *ptr;
            if ((t & 0x000FFFFFFFFFFFFF) && !(t & 0x7FF0000000000000)) *ptr = 0;
            ptr++;
        }
    }

    void audio_math::add_offset(float* p_buffer, float p_delta, size_t p_count) {
        for (size_t n = 0; n < p_count; ++n) p_buffer[n] += p_delta;
    }
    void audio_math::add_offset(double* p_buffer, double p_delta, size_t p_count) {
        for (size_t n = 0; n < p_count; ++n) p_buffer[n] += p_delta;
    }

	void audio_math::convert(const float* in, float* out, size_t count) {
        memcpy(out, in, count * sizeof(float));
    }
    void audio_math::convert(const float* in, float* out, size_t count, float scale) {
        audio_math::scale(in, count, out, scale);
    }
    void audio_math::convert(const double* in, double* out, size_t count) {
        memcpy(out, in, count * sizeof(double));
    }
    void audio_math::convert(const double* in, double* out, size_t count, double scale) {
        audio_math::scale(in, count, out, scale);
    }

    void audio_math::convert(const float* in, double* out, size_t count) {
        // optimize me
        noopt_convert(in, out, count);
    }
    void audio_math::convert(const float* in, double* out, size_t count, double scale) {
        // optimize me
        noopt_scale(in, count, out, scale);
    }
    void audio_math::convert(const double* in, float* out, size_t count) {
        // optimize me
        noopt_convert(in, out, count);
    }
    void audio_math::convert(const double* in, float* out, size_t count, double scale) {
        // optimize me
        noopt_scale(in, count, out, scale);
    }


    typedef char store24_t;
    static store24_t* store24(store24_t* out, int32_t in) {
        *(out++) = ((store24_t*)&in)[0];
        *(out++) = ((store24_t*)&in)[1];
        *(out++) = ((store24_t*)&in)[2];
        return out;
    }
    static store24_t* store24p(store24_t* out, int32_t in) {
        *(int32_t*)out = in;
        return out + 3;
    }

    static constexpr int32_t INT24_MAX = 0x7FFFFF, INT24_MIN = -0x800000;

    template<typename float_t> void convert_to_int24_noopt(float_t const* in, size_t count, void* out, float_t scale) {
        if (count == 0) return;
        --count;
        auto ptr = reinterpret_cast<store24_t*>(out);
        constexpr float_t lo = INT24_MIN, hi = INT24_MAX;
        while (count) {
            auto vf = *in++ * scale;
            if (vf < lo) vf = lo;
            else if (vf > hi) vf = hi;
            ptr = store24p(ptr, audio_math::rint32(vf));
            --count;
        }

        auto vf = *in * scale;
        if (vf < lo) vf = lo;
        else if (vf > hi) vf = hi;
        store24(ptr, audio_math::rint32(vf));
    }
#ifdef AUDIO_MATH_SSE
#if allowAVX
    static void f64_to_i24_avx(double const* in, size_t n, uint8_t* out, double scale) {
        const __m128i pi0 = _mm_set_epi8(-128, -128, -128, -128, 14, 13, 12, 10, 9, 8, 6, 5, 4, 2, 1, 0);
        const __m128i pi1 = _mm_set_epi8(4, 2, 1, 0, -128, -128, -128, -128, 14, 13, 12, 10, 9, 8, 6, 5);
        const __m128i pi2 = _mm_set_epi8(9, 8, 6, 5, 4, 2, 1, 0, -128, -128, -128, -128, 14, 13, 12, 10);
        const __m128i pi3 = _mm_set_epi8(14, 13, 12, 10, 9, 8, 6, 5, 4, 2, 1, 0, -128, -128, -128, -128);
        const auto mul = _mm256_set1_pd(scale);

        // PROBLEM: if we want to handle wildly out-of-bounds values, we can't do int clipping!
        // float clipping is sadly considerably slower than int clipping
        const auto lo = _mm256_set1_pd(INT24_MIN);
        const auto hi = _mm256_set1_pd(INT24_MAX);

        while (n >= 4 * 4) {
            auto f0 = _mm256_mul_pd(_mm256_loadu_pd(in + 0), mul);
            auto f1 = _mm256_mul_pd(_mm256_loadu_pd(in + 4), mul);
            auto f2 = _mm256_mul_pd(_mm256_loadu_pd(in + 8), mul);
            auto f3 = _mm256_mul_pd(_mm256_loadu_pd(in + 12), mul);
            f0 = _mm256_max_pd(_mm256_min_pd(f0, hi), lo);
            f1 = _mm256_max_pd(_mm256_min_pd(f1, hi), lo);
            f2 = _mm256_max_pd(_mm256_min_pd(f2, hi), lo);
            f3 = _mm256_max_pd(_mm256_min_pd(f3, hi), lo);
            __m128i w0 = _mm256_cvtpd_epi32(f0);
            __m128i w1 = _mm256_cvtpd_epi32(f1);
            __m128i w2 = _mm256_cvtpd_epi32(f2);
            __m128i w3 = _mm256_cvtpd_epi32(f3);

            // _mm_shuffle_epi8 : SSSE3
            w0 = _mm_shuffle_epi8(w0, pi0);
            w1 = _mm_shuffle_epi8(w1, pi1);
            w2 = _mm_shuffle_epi8(w2, pi2);
            w3 = _mm_shuffle_epi8(w3, pi3);

            // _mm_blend_epi16 : SSE4.1
            __m128i u0 = _mm_blend_epi16(w0, w1, 0xC0);
            __m128i u1 = _mm_blend_epi16(w1, w2, 0xF0);
            __m128i u2 = _mm_blend_epi16(w2, w3, 0xFC);

            _mm_storeu_si128((__m128i*)(out + 0), u0);
            _mm_storeu_si128((__m128i*)(out + 16), u1);
            _mm_storeu_si128((__m128i*)(out + 32), u2);

            in += 4 * 4;
            out += 16 * 3;
            n -= 4 * 4;
        }

        convert_to_int24_noopt(in, n, out, scale);
    }
#endif // allowAVX
    static void f64_to_i24_sse41(double const* in, size_t n, uint8_t* out, double scale) {
        const __m128i pi0 = _mm_set_epi8(-128, -128, -128, -128, 14, 13, 12, 10, 9, 8, 6, 5, 4, 2, 1, 0);
        const __m128i pi1 = _mm_set_epi8(4, 2, 1, 0, -128, -128, -128, -128, 14, 13, 12, 10, 9, 8, 6, 5);
        const __m128i pi2 = _mm_set_epi8(9, 8, 6, 5, 4, 2, 1, 0, -128, -128, -128, -128, 14, 13, 12, 10);
        const __m128i pi3 = _mm_set_epi8(14, 13, 12, 10, 9, 8, 6, 5, 4, 2, 1, 0, -128, -128, -128, -128);
        const auto mul = _mm_set1_pd(scale);

        // PROBLEM: if we want to handle wildly out-of-bounds values, we can't do int clipping!
        // float clipping is sadly considerably slower than int clipping
        const auto lo = _mm_set1_pd(INT24_MIN);
        const auto hi = _mm_set1_pd(INT24_MAX);

        while (n >= 4 * 4) {
            auto f0 = _mm_mul_pd(_mm_loadu_pd(in + 0), mul);
            auto f1 = _mm_mul_pd(_mm_loadu_pd(in + 2), mul);
            auto f2 = _mm_mul_pd(_mm_loadu_pd(in + 4), mul);
            auto f3 = _mm_mul_pd(_mm_loadu_pd(in + 6), mul);
            auto f4 = _mm_mul_pd(_mm_loadu_pd(in + 8), mul);
            auto f5 = _mm_mul_pd(_mm_loadu_pd(in + 10), mul);
            auto f6 = _mm_mul_pd(_mm_loadu_pd(in + 12), mul);
            auto f7 = _mm_mul_pd(_mm_loadu_pd(in + 14), mul);
            f0 = _mm_max_pd(_mm_min_pd(f0, hi), lo);
            f1 = _mm_max_pd(_mm_min_pd(f1, hi), lo);
            f2 = _mm_max_pd(_mm_min_pd(f2, hi), lo);
            f3 = _mm_max_pd(_mm_min_pd(f3, hi), lo);
            f4 = _mm_max_pd(_mm_min_pd(f4, hi), lo);
            f5 = _mm_max_pd(_mm_min_pd(f5, hi), lo);
            f6 = _mm_max_pd(_mm_min_pd(f6, hi), lo);
            f7 = _mm_max_pd(_mm_min_pd(f7, hi), lo);

            

            __m128i w0 = _mm_unpacklo_epi64(_mm_cvtpd_epi32(f0), _mm_cvtpd_epi32(f1));
            __m128i w1 = _mm_unpacklo_epi64(_mm_cvtpd_epi32(f2), _mm_cvtpd_epi32(f3));
            __m128i w2 = _mm_unpacklo_epi64(_mm_cvtpd_epi32(f4), _mm_cvtpd_epi32(f5));
            __m128i w3 = _mm_unpacklo_epi64(_mm_cvtpd_epi32(f6), _mm_cvtpd_epi32(f7));

            // _mm_shuffle_epi8 : SSSE3
            w0 = _mm_shuffle_epi8(w0, pi0);
            w1 = _mm_shuffle_epi8(w1, pi1);
            w2 = _mm_shuffle_epi8(w2, pi2);
            w3 = _mm_shuffle_epi8(w3, pi3);

            // _mm_blend_epi16 : SSE4.1
            __m128i u0 = _mm_blend_epi16(w0, w1, 0xC0);
            __m128i u1 = _mm_blend_epi16(w1, w2, 0xF0);
            __m128i u2 = _mm_blend_epi16(w2, w3, 0xFC);

            _mm_storeu_si128((__m128i*)(out + 0), u0);
            _mm_storeu_si128((__m128i*)(out + 16), u1);
            _mm_storeu_si128((__m128i*)(out + 32), u2);

            in += 4 * 4;
            out += 16 * 3;
            n -= 4 * 4;
        }

        convert_to_int24_noopt(in, n, out, scale);
    }
    static void f32_to_i24_sse41(float const* in, size_t n, uint8_t* out, float scale) {
        const __m128i pi0 = _mm_set_epi8(-128, -128, -128, -128, 14, 13, 12, 10, 9, 8, 6, 5, 4, 2, 1, 0);
        const __m128i pi1 = _mm_set_epi8(4, 2, 1, 0, -128, -128, -128, -128, 14, 13, 12, 10, 9, 8, 6, 5);
        const __m128i pi2 = _mm_set_epi8(9, 8, 6, 5, 4, 2, 1, 0, -128, -128, -128, -128, 14, 13, 12, 10);
        const __m128i pi3 = _mm_set_epi8(14, 13, 12, 10, 9, 8, 6, 5, 4, 2, 1, 0, -128, -128, -128, -128);
        const __m128 mul = _mm_set1_ps(scale);

        // PROBLEM: if we want to handle wildly out-of-bounds values, we can't do int clipping!
        // float clipping is sadly considerably slower than int clipping
        const auto lo = _mm_set1_ps(INT24_MIN);
        const auto hi = _mm_set1_ps(INT24_MAX);

        while (n >= 4 * 4) {
            auto f0 = _mm_mul_ps(_mm_loadu_ps(in + 0), mul);
            auto f1 = _mm_mul_ps(_mm_loadu_ps(in + 4), mul);
            auto f2 = _mm_mul_ps(_mm_loadu_ps(in + 8), mul);
            auto f3 = _mm_mul_ps(_mm_loadu_ps(in + 12), mul);
            f0 = _mm_min_ps(_mm_max_ps(f0, lo), hi);
            f1 = _mm_min_ps(_mm_max_ps(f1, lo), hi);
            f2 = _mm_min_ps(_mm_max_ps(f2, lo), hi);
            f3 = _mm_min_ps(_mm_max_ps(f3, lo), hi);
            __m128i w0 = _mm_cvtps_epi32(f0);
            __m128i w1 = _mm_cvtps_epi32(f1);
            __m128i w2 = _mm_cvtps_epi32(f2);
            __m128i w3 = _mm_cvtps_epi32(f3);

            // _mm_shuffle_epi8 : SSSE3
            w0 = _mm_shuffle_epi8(w0, pi0);
            w1 = _mm_shuffle_epi8(w1, pi1);
            w2 = _mm_shuffle_epi8(w2, pi2);
            w3 = _mm_shuffle_epi8(w3, pi3);

            // _mm_blend_epi16 : SSE4.1
            __m128i u0 = _mm_blend_epi16(w0, w1, 0xC0);
            __m128i u1 = _mm_blend_epi16(w1, w2, 0xF0);
            __m128i u2 = _mm_blend_epi16(w2, w3, 0xFC);

            _mm_storeu_si128((__m128i*)(out + 0), u0);
            _mm_storeu_si128((__m128i*)(out + 16), u1);
            _mm_storeu_si128((__m128i*)(out + 32), u2);

            in += 4 * 4;
            out += 16 * 3;
            n -= 4 * 4;
        }
        
        convert_to_int24_noopt(in, n, out, scale);
    }

#endif // AUDIO_MATH_SSE

    void audio_math::convert_to_int24(const float* in, size_t count, void* out, float scale) {
        scale *= 0x800000;
        
#ifdef AUDIO_MATH_SSE
        if (haveSSE41) {
            f32_to_i24_sse41(in, count, (uint8_t*)out, scale); return;
        }
#endif
        convert_to_int24_noopt(in, count, out, scale);
    }
    void audio_math::convert_to_int24(const double* in, size_t count, void* out, double scale) {
        scale *= 0x800000;
#ifdef AUDIO_MATH_SSE
#if allowAVX
        if (haveAVX) {
            f64_to_i24_avx(in, count, (uint8_t*)out, scale); return;
        }
#endif // allowAVX
        if (haveSSE41) {
            f64_to_i24_sse41(in, count, (uint8_t*)out, scale); return;
        }
#endif // AUDIO_MATH_SSE
        convert_to_int24_noopt(in, count, out, scale);
    }

}