#include <iostream>
#include <immintrin.h>
// author: Jan Lemeire
// This code shows how to reduce an array with vector instructions
//
// OUTCOMMENT THE cout's and preceding stores TO USE IT!!!!!!!!!!!!!!!!!

using namespace std;

float sumOf16Floats(float indata[16]){
	__m256 ymm0, ymm1, ymm2, ymm3, ymm4;
	float outdata[8];

	ymm0 = _mm256_loadu_ps(indata);
	ymm1 = _mm256_loadu_ps(indata + 8);
	ymm2 = _mm256_hadd_ps(ymm0, ymm1); // see https://www.officedaytime.com/simd512e/simdimg/si.php?f=haddps

	_mm256_storeu_ps(outdata, ymm2);
	cout << "Out vector : " << outdata[0] << " " << outdata[1] << " " << outdata[2] << " " << outdata[3] << " " << outdata[4] << " " << outdata[5] << " " << outdata[6] << " " << outdata[7] << endl;

	ymm3 = _mm256_hadd_ps(ymm2, ymm2);

	_mm256_storeu_ps(outdata, ymm3);
	cout << "Out vector2: " << outdata[0] << " " << outdata[1] << " " << outdata[2] << " " << outdata[3] << " " << outdata[4] << " " << outdata[5] << " " << outdata[6] << " " << outdata[7] << endl;

	ymm4 = _mm256_hadd_ps(ymm3, ymm3);
	_mm256_storeu_ps(outdata, ymm4);

	cout << "Out vector3: " << outdata[0] << " " << outdata[1] << " " << outdata[2] << " " << outdata[3] << " " << outdata[4] << " " << outdata[5] << " " << outdata[6] << " " << outdata[7] << endl;

	float sum_v = outdata[0] + outdata[4]; // last two sub results 
	cout << "Sum intrinsics = " << sum_v << endl;

	float sum = 0;
	for (int i = 0; i < 16; i++)
		sum += indata[i];
	cout << "Sum pure = " << sum << endl;
	return sum_v;
}

double sumOf8Doubles(double indata[8]) {
	__m256d ymm0, ymm1, ymm2, ymm3, ymm4;
	double outdata[4];

	ymm0 = _mm256_loadu_pd(indata);
	ymm1 = _mm256_loadu_pd(indata+4);
	ymm2 = _mm256_hadd_pd(ymm0, ymm1); // https://www.officedaytime.com/simd512e/simdimg/si.php?f=haddps

	_mm256_storeu_pd(outdata, ymm2);
	cout << "Out vector: " << outdata[0] << " " << outdata[1] << " " << outdata[2] << " " << outdata[3] << endl;

	ymm3 = _mm256_hadd_pd(ymm2, ymm2);

	_mm256_storeu_pd(outdata, ymm3);
	cout << "Out vector: " << outdata[0] << " " << outdata[1] << " " << outdata[2] << " " << outdata[3] << endl;

	double sum_d = outdata[0] + outdata[2];   // last two sub results 
	
	cout << "sum vectorized: " << sum_d << endl;

	//ymm3 = _mm256_permute2f128_pd(ymm2, ymm2, 0x01); // swap again
	//ymm4 = _mm256_add_pd(ymm2, ymm3);
	//_mm256_storeu_pd(outdata, ymm4);
	//cout << "Out vector: " << outdata[0] << " " << outdata[1] << " " <<outdata[2] << " " << outdata[3] << endl;

	double sum = 0;
	for (int i = 0; i < 8; i++)
		sum += indata[i];
	cout << "Sum pure = " << sum << endl;
	return sum_d;
}

int main()
{
	
	float indata[16] = { 1.0, 2.0, 3.0, 5.0, 5.0, 6.0, 7.0, 18.0, 100, 101, 102, 105, 200, 201, 202, 209 };
	sumOf16Floats(indata);

	double indata2[8] = { 1.0, 3.0, 4.0, 5.0 , 12.0, 13.0, 12.0, 15.0 };
	sumOf8Doubles(indata2);
	

	cout << endl << "Press ENTER to close window...";
	char c = cin.get();
}

// the following code does a sum of 4 doubles, in a rather strange way... and I think it's not faster than just scalar
//   taken from here: https://software.intel.com/en-us/blogs/2016/04/12/can-you-write-a-vectorized-reduction-operation

// sum of doubles
//double indata[4] = { 2.0, 3.0f, 2.0f, 5.0f };
//double outdata[4] = { 0.0, 0.0, 0.0, 0.0 };
//
//__m256d ymm0, ymm1, ymm2, ymm3, ymm4;
//
//int main()
//{
//	ymm0 = _mm256_loadu_pd(indata);
//	ymm1 = _mm256_permute_pd(ymm0, 0x05); // first two elements are swapped, second two as well
//	ymm2 = _mm256_add_pd(ymm0, ymm1);
//	ymm3 = _mm256_permute2f128_pd(ymm2, ymm2, 0x01); // swap again
//	ymm4 = _mm256_add_pd(ymm2, ymm3);
//	_mm256_storeu_pd(outdata, ymm4);
//	cout << "Out vector: " << outdata[0] << " " << outdata[1] << " " <<outdata[2] << " " << outdata[3] << endl;
//
//	cout << endl << "Press ENTER to close window...";
//	char c = cin.get();
//}