#include <iostream>

#include <x86intrin.h>

// Return horizontal sum for 128 bit packed float vector
float hsum_ps_sse3(__m128 val) 
{
    __m128 shuf = _mm_movehdup_ps(val); // broadcast elements 3,1 to 2,0
    __m128 sums = _mm_add_ps(val, shuf);
    shuf        = _mm_movehl_ps(shuf, sums);    // high half -> low half
    sums        = _mm_add_ss(sums, shuf);
    return        _mm_cvtss_f32(sums);
}

// Return horizontal sum for 256 bit packed float vector
// https://stackoverflow.com/a/35270026/73878 -- Peter Cordes is da man
float hsum256_ps_avx(__m256 val) 
{
    __m128 vlow  = _mm256_castps256_ps128(val);
    __m128 vhigh = _mm256_extractf128_ps(val, 1);   // high 128
           vlow  = _mm_add_ps(vlow, vhigh);         // add the low 128
    return hsum_ps_sse3(vlow);  // and inline the sse3 version, which is optimal for AVX
}

void print_debug(__m256 val)
{
    float temp[8];
    
    _mm256_store_ps(temp, val);
    std::cout << "DEBUG: " 
              << temp[0] << " -- " << temp[1] << " -- "
              << temp[2] << " -- " << temp[3] << " -- "
              << temp[4] << " -- " << temp[5] << " -- "
              << temp[6] << " -- " << temp[7] << std::endl;
}

// Return max value across (i.e. horizontally) vector
// Algorithm came from here: https://stackoverflow.com/a/9798369/73878
float max_ps_avx(__m256 val)
{
    __m256 m0 = _mm256_permute2f128_ps(val, val, 1); // Permute 128-bit values   
    __m256 m1 = _mm256_max_ps(val, m0); // Compute max of 4 sets
    
    __m256 m2 = _mm256_permute_ps(m1, 0x4E); // Permute max values
    __m256 m3 = _mm256_max_ps(m1, m2); // Compute max of 2 sets
    
    __m256 m4 = _mm256_permute_ps(m3, 0x3B); // Permute max values    
    __m256 m  = _mm256_max_ps(m3, m4);  // Compute (final) max
    
    return _mm_cvtss_f32(_mm256_extractf128_ps(m, 0));
}

// Return min value across (i.e. horizontally) vector
// Algorithm came from here: https://stackoverflow.com/a/9798369/73878
float min_ps_avx(__m256 val)
{
    __m256 m0 = _mm256_permute2f128_ps(val, val, 1); // Permute 128-bit values   
    __m256 m1 = _mm256_min_ps(val, m0); // Compute min of 4 sets
    
    __m256 m2 = _mm256_permute_ps(m1, 0x4E); // Permute max values
    __m256 m3 = _mm256_min_ps(m1, m2); // Compute min of 2 sets
    
    __m256 m4 = _mm256_permute_ps(m3, 0x3B); // Permute max values    
    __m256 m  = _mm256_min_ps(m3, m4);  // Compute (final) min
    
    return _mm_cvtss_f32(_mm256_extractf128_ps(m, 0));
}

// g++ -std=c++11 -O3 -mavx -o ps_test ps_test.cc
int main(int argc, char** argv)
{
    std::cout << "Packed Float (PS) Test\n" << std::endl;
    
    float raw_data[8] = {1.0, 2.0, -3.0, 4.0, 5.0, -6.0, 7.0 , 8.0};
    
    // Generate signmask
    __m256 SIGNMASK = _mm256_castsi256_ps(_mm256_set1_epi32(0x80000000));
    
    // Load data
    __m256 data = _mm256_loadu_ps(raw_data);
    
    // Take absolute value of data
    __m256 abs_data = _mm256_andnot_ps(SIGNMASK, data);
    
    // Sum data
    float total = hsum256_ps_avx(data);
    float total_abs = hsum256_ps_avx(abs_data);
    
    std::cout << "Total: " << total << " -- " << total_abs << std::endl;
    std::cout << "  Max: " << max_ps_avx(abs_data) << std::endl;
    std::cout << "  Min: " << min_ps_avx(abs_data) << std::endl;
    
    return 0;
}