#include <stdio.h>
#include <stdlib.h>

struct vector
{
    float data[4];
};

float sumVec(struct vector v){
    return v.data[0]+v.data[1]+v.data[2]+v.data[3];
}

float det(struct vector matrix[4])
{
    // matrix - tablica czterech wektorów w komentarzach kolejno A, B, C i D
    struct vector m;
    float result;

    /* Obliczenia do wykonania:
    det =   (a0*b1-a1*b0)*(c2*d3-c3*d2) +
        (a2*b0-a0*b2)*(c1*d3-c3*d1) +
        (a0*b3-a3*b0)*(c1*d2-c2*d1) +
        (a1*b2-a2*b1)*(c0*d3-c3*d0) +

        (a2*b3-a3*b2)*(c0*d1-c1*d0) +
        (a3*b1-a1*b3)*(c0*d2-c2*d0)

    Jak widać w druga suma jest zbyt mała by uzupełnić nią 4-elementowy wektor
    Stosujemy zatem następujące przekształcenia:

    det =   (a0*b1-a1*b0)*(c2*d3-c3*d2) +
        (a2*b0-a0*b2)*(c1*d3-c3*d1) +
        (a0*b3-a3*b0)*(c1*d2-c2*d1) +
        (a1*b2-a2*b1)*(c0*d3-c3*d0) +

        (a2*b3-a3*b2)*(c0*d1-c1*d0) +
        (a3*b1-a1*b3)*(c0*d2-c2*d0) +
              0      *      0       +
              0      *      0

    det =   (a0*b1-a1*b0)*(c2*d3-c3*d2) +
        (a2*b0-a0*b2)*(c1*d3-c3*d1) +
        (a0*b3-a3*b0)*(c1*d2-c2*d1) +
        (a1*b2-a2*b1)*(c0*d3-c3*d0) +

        (a2*b3-a3*b2)*(c0*d1-c1*d0) +
        (a3*b1-a1*b3)*(c0*d2-c2*d0) +
        (a0*b0-a0*b0)*(c0*d0-c0*d0) +
        (a0*b0-a0*b0)*(c0*d0-c0*d0)

    Teraz kolumny indeksów elementów wektora mogą stać się argumentami kolejnych wywołań shufps.
    Każdy wektor przestawimy zgodnie z indeksami a następnie będziemy wykonywać operacje
    na wektorach czteroelementowych.
    */

    asm (
        /* Obliczenie 1 składowego wyniku
        (a0*b1-a1*b0)
        (a2*b0-a0*b2)
        (a0*b3-a3*b0)
        (a1*b2-a2*b1)
         A' B' A" B"
        */
        "movups (%0), %%xmm0     \n\t"
        "movups (%1), %%xmm1 \n\t"

        "shufps $0x48 , %%xmm0, %%xmm0 \n\t"            //A' (1, 0, 2, 0)
        "shufps $0xb1 , %%xmm1, %%xmm1 \n\t"            //B' (2, 3, 0, 1)
        "mulps %%xmm0,%%xmm1 \n\t"              //A'*B'

        // ponowne ustawienie zawartości rejestrów
        // Obliczenie iloczynów z prawej strony
        "movups (%0), %%xmm0 \n\t"
        "movups (%1), %%xmm2 \n\t"

        "shufps $0xb1,%%xmm0,%%xmm0 \n\t"           //A" (2, 3, 0, 1)
        "shufps $0x48,%%xmm2,%%xmm2 \n\t"           //B" (1, 0, 2, 0)
        "mulps  %%xmm0,%%xmm2 \n\t"                 //A"*B"

        // Wyliczenie różnicy
        "subps %%xmm1, %%xmm2 \n\t"                 //A'*B'-A"*B" po wektoryzacji

        /*
        (c2*d3-c3*d2)
        (c1*d3-c3*d1)
        (c1*d2-c2*d1)
        (c0*d3-c3*d0)
         C' D' C" D"
        */
        "movups (%2), %%xmm0     \n\t"              //C,
        "movups (%3), %%xmm1 \n\t"                  //D

        "shufps $0x16 , %%xmm0, %%xmm0 \n\t"            //C' (0, 1, 1, 2)
        "shufps $0xef , %%xmm1, %%xmm1 \n\t"            //D' (3, 2, 3, 3)
        "mulps %%xmm0,%%xmm1 \n\t"              // C'*D'

        "movups (%2), %%xmm0 \n\t"              //C
        "movups (%3), %%xmm3 \n\t"                  //D

        "shufps $0xef,%%xmm0 ,%%xmm0 \n\t"              //C" (3, 2, 3, 3)
        "shufps $0x16, %%xmm3, %%xmm3 \n\t"             //D" (0, 1, 1, 2)
        "mulps  %%xmm0,%%xmm3 \n\t"                 //C"*D"

        "subps %%xmm1, %%xmm3 \n\t"                 //C'*D'-C"*D"

        "mulps %%xmm2,%%xmm3 \n\t"              //(A'*B'-A"*B")*(C'*D'-C"*D")  - wynik skladowy 1


        /* Obliczenie 2 składowego wyniku
        (a2*b3-a3*b2)
        (a3*b1-a1*b3)
        (a0*b0-a0*b0) <- 0
        (a0*b0-a0*b0) <- 0
         A' B' A" B"
        */
        "movups (%0), %%xmm0     \n\t"              //A
        "shufps $0xe , %%xmm0, %%xmm0  \n\t"            //A' (0, 0, 3, 2)

        "movups (%1), %%xmm2 \n\t"                  //B
        "shufps $0x7 , %%xmm2, %%xmm2 \n\t"             //B' (0, 0, 1, 3)

        "mulps %%xmm0, %%xmm2 \n\t"                 //A'*B'

        "movups (%0), %%xmm0     \n\t"              //A
        "shufps $0x7 , %%xmm0, %%xmm0 \n\t"             //A" (0, 0, 1, 3)

        "movups (%1), %%xmm4     \n\t"              //B
        "shufps $0xe , %%xmm4, %%xmm4 \n\t"             //B" (0, 0, 3, 2)

        "mulps %%xmm0,%%xmm4 \n\t"              //A"*B"

        "subps %%xmm4, %%xmm2 \n\t"                 //A'*B'-A"*B"

        /*
        (c0*d1-c1*d0)
        (c0*d2-c2*d0)
        (c0*d0-c0*d0) <- 0
        (c0*d0-c0*d0) <- 0
         C' D' C" D"
        */
        "movups (%2), %%xmm0 \n\t"              //C
        "shufps $0x0 , %%xmm0, %%xmm0 \n\t"         //C' (0, 0, 0, 0)

        "movups (%3), %%xmm4 \n\t"                  //D
        "shufps $0x9 , %%xmm4, %%xmm4 \n\t"             //D' (0, 0, 2, 1)

        "mulps %%xmm0, %%xmm4 \n\t"                 //C * D

        "movups (%2), %%xmm0 \n\t"              //C
        "shufps $0x9 , %%xmm0, %%xmm0 \n\t"         //C" (0, 0, 2, 1)

        "movups (%3), %%xmm5 \n\t"                  //D
        "shufps $0x0 , %%xmm5, %%xmm5 \n\t"         //D" (0, 0, 0, 0)

        "mulps %%xmm5,%%xmm0 \n\t"              //C" * D"

        "subps %%xmm4, %%xmm0 \n\t"                 //C'*D'- C"*D"

        "mulps %%xmm0,%%xmm2 \n\t"              //skladowa wyniku 2

        "subps %%xmm2,%%xmm3 \n\t"              //roznica skladowych = wynik
        "movups %%xmm3,(%4) \n\t"               //wektor rozwiazan
        :
        : "r" (matrix[0].data),"r" (matrix[1].data) , "r" (matrix[2].data), "r"(matrix[3].data), "r" (m.data)
    );

    int i;
    printf("Sumy częściowe: \n");
    for(i = 0; i<4; i++)
        printf("%f \n",m.data[i]);

    return sumVec(m);
};


int main()
{

    float deter;

    struct vector matrix[4];

    // tablica użyta jedynie do ułatwienia zapisu macierzowego
    // matrix można uzupełniać przez matrix[i].data[j] bez inicjalizacji tablicy
    float m[4][4] = {{1, 2, 3, 4},
             {5, 6, 7, 8},
             {9, 1, 2, 3},
             {4, 5, 6, 8} };

    int i, j;

    for(i = 0; i < 4; i++){
        for(j = 0; j < 4;  j++)
            matrix[i].data[j] = m[i][j];
    }


    deter=det(matrix);
    printf("Wyznacznik: \n");
    printf("%f \n",deter);

    return 0;
}