// Compilation (linux) : g++ CrossProd.cpp -O3 -o prog

#include <iostream>
#include <boost/timer.hpp>
#include <iostream>

using namespace std;

const int N = 10;
const int n = 10000;
const int p = 200;

double loop4()
{
    double **x = new double*[n];
    for (int i = 0; i < n; i++)
    {
        x[i] = new double[p];
        for (int j = 0; j < p ; j++)
            x[i][j] = (i*j) % 2;
    }
        
    double **mat = new double*[p];
    for (int j = 0; j < p; j++)
    {
        mat[j] = new double[p];
        std::fill(mat[j], mat[j] + p, 0.0);
    }
        
    int nbOp = 0;

    boost::timer t;
    for (int i = 0; i < n ; i++)
    {
        for (int j1 = 0; j1 < p ; j1++)
        {
            double xij1 = x[i][j1];
            for (int j2 = 0; j2 < p ; j2+=16)
            {
                for (int k = 0 ; (k < 16) && (j2+k < p); k++)
                {
                    mat[j1][j2+k] = mat[j1][j2+k] + xij1 * x[i][j2+k];
                }
                nbOp = nbOp + 16;
            }
        }
    }
    double d = t.elapsed();
    
    for (int j = 0; j < p ; j++)    
        delete [] mat[j];
    delete [] mat;

    for (int i = 0; i < n ; i++) 
        delete [] x[i];
    delete [] x;
    
    return d;
}

double loop3()
{
    double **x = new double*[n];
    for (int i = 0; i < n; i++)
    {
        x[i] = new double[p];
        for (int j = 0; j < p ; j++)
            x[i][j] = (i*j) % 2;
    }
        
    double **mat = new double*[p];
    for (int j = 0; j < p; j++)
    {
        mat[j] = new double[p];
        std::fill(mat[j], mat[j] + p, 0.0);
    }
    
    int nbOp = 0;

    boost::timer t;
    for (int i = 0; i < n ; i++)
    {
        for (int j1 = 0; j1 < p ; j1++)
        {
            double xij1 = x[i][j1];
            for (int j2 = 0; j2 < p ; j2++)
            {
                mat[j1][j2] = mat[j1][j2] + xij1 * x[i][j2];
                nbOp = nbOp + 1;
            }
        }
    }
    double d = t.elapsed();
    
    for (int j = 0; j < p ; j++)    
        delete [] mat[j];
    delete [] mat;

    for (int i = 0; i < n ; i++) 
        delete [] x[i];
    delete [] x;
    
    return d;
}

double loop2()
{
    double **x = new double*[n];
    for (int i = 0; i < n; i++)
    {
        x[i] = new double[p];
        for (int j = 0; j < p ; j++)
            x[i][j] = (i*j) % 2;
    }
        
    double **mat = new double*[p];
    for (int j = 0; j < p; j++)
    {
        mat[j] = new double[p];
        std::fill(mat[j], mat[j] + p, 0.0);
    }
        
        
    int nbOp = 0;

    boost::timer t;
    for (int i = 0; i < n ; i++)
    {
        for (int j1 = 0; j1 < p ; j1++)
        {
            for (int j2 = 0; j2 < p ; j2++)
            {
                mat[j1][j2] = mat[j1][j2] + x[i][j1] * x[i][j2];
                nbOp = nbOp + 1;
            }
        }
    }
    double d = t.elapsed();
    
    for (int j = 0; j < p ; j++)    
        delete [] mat[j];
    delete [] mat;

    for (int i = 0; i < n ; i++) 
        delete [] x[i];
    delete [] x;
    
    return d;
}

double loop1()
{
    double **x = new double*[n];
    for (int i = 0; i < n; i++)
    {
        x[i] = new double[p];
        for (int j = 0; j < p ; j++)
            x[i][j] = (i*j) % 2;
    }

    double **mat = new double*[p];
    for (int j = 0; j < p; j++) 
        mat[j] = new double[p];
        
    int nbOp = 0;

    boost::timer t;    
    for (int i = 0; i < n ; i++)
    {
        for (int j1 = 0; j1 < p ; j1++)
        {
            for (int j2 = 0; j2 < p ; j2++)
            {
                mat[j1][j2] = mat[j1][j2] + x[i][j1] * x[i][j2];
                nbOp = nbOp + 1;
            }
        }
    }
    double d = t.elapsed();
        
    for (int j = 0; j < p ; j++)    
        delete [] mat[j];
    delete [] mat;

    for (int i = 0; i < n ; i++) 
        delete [] x[i];
    delete [] x;
    
    return d;
}

double loop()
{
    double **x = new double*[n];
    for (int i = 0; i < n; i++)
    {
        x[i] = new double[p];
        for (int j = 0; j < p ; j++)
            x[i][j] = (i*j) % 2;
    }

    double **mat = new double*[p];
    for (int j = 0; j < p; j++) 
        mat[j] = new double[p];
        
    int nbOp = 0;

    boost::timer t;
    for (int j1 = 0; j1 < p ; j1++)
    {
        for (int j2 = 0; j2 < p ; j2++)
        {
            for (int i = 0; i < n ; i++)
            {
                mat[j1][j2] = mat[j1][j2] + x[i][j1] * x[i][j2];
                nbOp = nbOp + 1;
            }
        }
    }
    double d = t.elapsed();
    
    for (int j = 0; j < p ; j++)    
        delete [] mat[j];
    delete [] mat;

    for (int i = 0; i < n ; i++) 
        delete [] x[i];
    delete [] x;
    
    return d;
}
    
int main(int argc, char* argv[])
{
    cout << "debut des calculs..." << endl;
    { double sum = 0.0; for (int i = 0; i < N; ++i) { sum += loop(); } ; std::cout << "version 0: " << sum/N << std::endl; }
    { double sum = 0.0; for (int i = 0; i < N; ++i) { sum += loop1(); } ; std::cout << "version 1: " << sum/N << std::endl; }
    { double sum = 0.0; for (int i = 0; i < N; ++i) { sum += loop2(); } ; std::cout << "version 2: " << sum/N << std::endl; }
    { double sum = 0.0; for (int i = 0; i < N; ++i) { sum += loop3(); } ; std::cout << "version 3: " << sum/N << std::endl; }
    { double sum = 0.0; for (int i = 0; i < N; ++i) { sum += loop4(); } ; std::cout << "version 4: " << sum/N << std::endl; }             
    return 0;

}