#include <iostream>
#include <array>
#include <amp.h>
#include <stdio.h>
#include <tchar.h>
#include <immintrin.h>
using namespace concurrency;
using namespace std;

union ieee754_QNAN
{
    const float f;
    struct
    {
        const unsigned int mantissa:23, exp:8, sign:1;
    };

    ieee754_QNAN() : f(0.0), mantissa(0x7FFFFF), exp(0xFF), sign(0x0) {}
};

_MM_ALIGN16 const ieee754_QNAN absMask;
static const __m256 abs4Mask = _mm256_set1_ps(absMask.f);

const __m256 _mm256_abs_ps(const __m256& x)
{
    return _mm256_and_ps(abs4Mask, x);
}

int _tmain(int argc, _TCHAR* argv[])
{
    LARGE_INTEGER timeStart;
    QueryPerformanceCounter(&timeStart);

    const int wantedNumPrimes = 10000;
    const int step = 10000;

    std::vector<int> primes;
    std::vector<float> primesInv;
    primes.reserve(wantedNumPrimes * 2);
    primesInv.reserve(wantedNumPrimes * 2);

    primes.push_back(2);
    primesInv.push_back(1.0f / 2.0f);

    primes.push_back(3);
    primesInv.push_back(1.0f / 3.0f);

    primes.push_back(5);
    primesInv.push_back(1.0f / 5.0f);

    primes.push_back(7);
    primesInv.push_back(1.0f / 7.0f);

    primes.push_back(11);
    primesInv.push_back(1.0f / 11.0f);

    primes.push_back(13);
    primesInv.push_back(1.0f / 13.0f);

    primes.push_back(17);
    primesInv.push_back(1.0f / 17.0f);

    primes.push_back(19);
    primesInv.push_back(1.0f / 19.0f);

    int blockStart = 20;
    while(primes.size() < wantedNumPrimes)
    {
        int from = blockStart;
        int to = blockStart + step;
        blockStart = to;

        for(int n = from; n < to; n += 8)
        {
            __m256 n8 = _mm256_set_ps(n, n + 1, n + 2, n + 3, n + 4, n + 5, n + 6, n + 7);

            __m256 isPrimeMask8 = _mm256_setzero_ps();
            int isPrimeMask = 0;
            int wall = n + 7;
            for(int i = 0; i < primes.size(); i++)
            {
                if(primes[i] * primes[i] >= wall)
                    break;

                __m256 rcp = _mm256_set1_ps(primesInv[i]);

                __m256 divisionRes = _mm256_mul_ps(n8, rcp);
                __m256i rounded = _mm256_cvtps_epi32(divisionRes);

                __m256 rounded8 = _mm256_cvtepi32_ps(rounded);
                __m256 diff = _mm256_abs_ps(_mm256_sub_ps(rounded8, divisionRes));
                __m256 cmpResult = _mm256_cmp_ps(diff, _mm256_set1_ps(0.001f), _CMP_LT_OS);
                isPrimeMask8 = _mm256_or_ps(isPrimeMask8, cmpResult);

                isPrimeMask = _mm256_movemask_ps(isPrimeMask8);

                if(isPrimeMask == 255)
                    break;
            }

            if(isPrimeMask == 255)
                continue;

            if((isPrimeMask & 0x80) == 0)
            {
                int prime = n;
                primes.push_back(prime);
                primesInv.push_back(1.0f / prime);
            }

            if((isPrimeMask & 0x40) == 0)
            {
                int prime = n + 1;
                primes.push_back(prime);
                primesInv.push_back(1.0f / prime);
            }

            if((isPrimeMask & 0x20) == 0)
            {
                int prime = n + 2;
                primes.push_back(prime);
                primesInv.push_back(1.0f / prime);
            }

            if((isPrimeMask & 0x10) == 0)
            {
                int prime = n + 3;
                primes.push_back(prime);
                primesInv.push_back(1.0f / prime);
            }

            if((isPrimeMask & 0x8) == 0)
            {
                int prime = n + 4;
                primes.push_back(prime);
                primesInv.push_back(1.0f / prime);
            }

            if((isPrimeMask & 0x4) == 0)
            {
                int prime = n + 5;
                primes.push_back(prime);
                primesInv.push_back(1.0f / prime);
            }

            if((isPrimeMask & 0x2) == 0)
            {
                int prime = n + 6;
                primes.push_back(prime);
                primesInv.push_back(1.0f / prime);
            }

            if((isPrimeMask & 0x1) == 0)
            {
                int prime = n + 7;
                primes.push_back(prime);
                primesInv.push_back(1.0f / prime);
            }
        }
    }

    LARGE_INTEGER timeEnd;
    QueryPerformanceCounter(&timeEnd);
    float timeDiff = timeEnd.LowPart - timeStart.LowPart;

    LARGE_INTEGER freqPerSecond;
    QueryPerformanceFrequency(&freqPerSecond);
    float time = (timeDiff * 1000.0f)/(float)freqPerSecond.LowPart;

    for(int i = 0; i < 10000; i++)
    {
        cout << primes[i] << " ";
    }

    cout<<"Time elapsed: "<<time<<"ms";

    int a;
    cin>>a;

    return 0;
}