import std.datetime, std.stdio, std.random;
import ldc.attributes : fastmath;

@fastmath void compute_neq(float[] neq,
                           const float[] ux,
                           const float[] uy,
                           const float[] rho,
                           const float[] ex,
                           const float[] ey,
                           const float[] w,
                           const size_t N)
{
    for (size_t idx = 0; idx < N * N; ++idx) {
        float usqr = ux[idx] * ux[idx] + uy[idx] * uy[idx];
        for (size_t q = 0; q < 9; ++q) {
            float eu = 3.0f * (ex[q] * ux[idx] + ey[q] * uy[idx]);
            float tmp = 1.0f + eu + 0.5f * eu * eu - 1.5f * usqr;
            tmp *= w[q] * rho[idx];
            neq[idx * 9 + q] = tmp;
        }
    }
}


void main(string[] args)
{
    immutable size_t N = 1000;
    immutable int length = N * N;
    immutable int size = float.sizeof * length;
    immutable size_t benchmarkCount = 1000;
    const float[9] ex = [0., 1., 0., -1., 0., 1., -1., -1., 1.];
    const float[9] ey = [0., 0., 1., 0., -1., 1., 1., -1., -1.];
    const float[9] w = [4./9., 1./9., 1./9., 1./9., 1./9., 1./36., 1./36., 1./36., 1./36.];
    float[] neq = new float[N*N*9];
    foreach(int i,e; neq) neq[i] = 0.0f;
    float[] ux = new float[N*N];
    float[] uy = new float[N*N];
    float[] rho = new float[N*N];
    Random gen;
    foreach(int i,e; ux) ux[i] = uniform(0.0f, 1.0f, gen);
    foreach(int i,e; uy) uy[i] = uniform(0.0f, 1.0f, gen);
    foreach(int i,e; rho) rho[i] = 1.0f;
    StopWatch sw;
    sw.start();
    for (uint c=0; c<benchmarkCount; ++c)
        compute_neq(neq, ux, uy, rho, ex, ey, w, N);
    sw.stop();
    float t = sw.peek().usecs;
    writef("neq[0]=%f\tt=%f ms", neq[0], t / 1000.0f / benchmarkCount);
}