#include <iostream>
#include <ctime>
#include <cmath>
#include <GL\freeglut.h>
#include <omp.h>
#include <amp.h>
#include <amp_math.h>
using namespace concurrency;

const int N = 15000;
int button = 2;

float POSX[N];
float POSY[N];
float POSR[N + N];
float AX_M[N];
float AY_M[N];


void initialization()
{
    glClearColor(0, 0, 0, 1);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluOrtho2D(-1000, 1000, -1000, 1000);
}


void display()
{
    glClear(GL_COLOR_BUFFER_BIT);
    glPushMatrix();
    glBegin(GL_POINTS);

    for (int draw = 0; draw < N; draw++) {
        glColor3f(1, 1, 1);
        glVertex2f(POSX[draw] += POSR[draw], POSY[draw] += POSR[draw + N]);
    }

    glEnd();

    glPopMatrix();
    glutSwapBuffers();
}


void OMP_Math()
{
    int den;
#pragma omp parallel for private(den, AX_M, AY_M)
    for (int first = 0; first < N; first++) {
        for (int next = 0; next < N; next++) {
            if ((first != next) && ((POSX[next] != POSX[first]) || (POSY[next] != POSY[first]))) {
                den = 10 * sqrt((POSX[next] - POSX[first])*(POSX[next] - POSX[first]) + (POSY[next] - POSY[first])*(POSY[next] - POSY[first]));
                AX_M[first] = (POSX[next] - POSX[first]) / den;
                AY_M[first] = (POSY[next] - POSY[first]) / den;
            }
        }
        POSR[first] += AX_M[first];
        POSR[first + N] += AY_M[first];
    }
}

void APM_Math()
{
    array_view<float, 1> Matrix_X(N, POSX), Matrix_Y(N, POSY), Matrix_R(N + N, POSR), AX(N, AX_M), AY(N, AY_M);
    Matrix_R.discard_data();
    parallel_for_each(
        Matrix_R.extent, [=](index<1> first) restrict(amp) {
        for (int next = 0; next < N; next++) {
            if ((first[0] != next) && ((Matrix_X[next] != Matrix_X[first[0]]) || (Matrix_Y[next] != Matrix_Y[first[0]]))) {
                int den;
                den = 10 * fast_math::sqrt((Matrix_X[next] - Matrix_X[first[0]])*(Matrix_X[next] - Matrix_X[first[0]]) + (Matrix_Y[next] - Matrix_Y[first[0]])*(Matrix_Y[next] - Matrix_Y[first[0]]));
                AX[first[0]] = (Matrix_X[next] - Matrix_X[first[0]]) / den;
                AY[first[0]] = (Matrix_Y[next] - Matrix_Y[first[0]]) / den;
            }
        }
        Matrix_R[first[0]] += AX[first[0]];
        Matrix_R[first[0] + N] += AY[first[0]];
    });
    Matrix_R.synchronize();
}

void Timer(int)
{
    if (button != 1)
    {
        OMP_Math();
    }

    else
    {
        APM_Math();
    }

    display();
    glutTimerFunc(0, Timer, 0);
}


void characters()
{
    for (int i = 0; i < N; i++) {
        POSX[i] = rand() % 1000 - 500;
        POSY[i] = rand() % 1000 - 500;

        for (int j = 0; j < i; j++) {
            if ((i != j) && (POSX[i] == POSX[j]) && (POSY[i] == POSY[j]))
                i--;
        }
    }
}


void keyboard(unsigned char key, int x, int y)
{
    switch (key)
    {
    case '1': button = 1;
        break;
    case '2': button = 2;
        break;
    }

}


int main(int argc, char **argv)
{
    srand(time(0));
    setlocale(LC_ALL, "russian");

    characters();

    std::cout << "[1] -- Select GPU\n";
    std::cout << "[2] -- Select CPU\n";

    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
    glutInitWindowPosition(100, 100);
    glutInitWindowSize(1000, 1000);
    glutCreateWindow("YOBA N-Body.Ver. 0.7a");
    glutDisplayFunc(display);

    glutTimerFunc(0, Timer, 0);
    glutKeyboardFunc(keyboard);
    initialization();
    glutMainLoop();

    return 0;
}