struct float4x4
{
    float4 cols[4];
};

constexpr float4x4 Identity4x4 = {
    float4{1, 0, 0, 0},
    float4{0, 1, 0, 0},
    float4{0, 0, 1, 0},
    float4{0, 0, 0, 1},
};

float4x4
operator*(const float4x4& a, const float4x4& b)
{
    float4x4 ret {};
    for (uint64_t row = 0; row < 4; ++row) {
        for (uint64_t col = 0; col < 4; ++col) {
            for (uint64_t k = 0; k < 4; ++k)
                ((float*)&ret.cols[col])[row] += ((float*)&a.cols[k])[row] * ((float*)&b.cols[col])[k];
        }
    }
    return ret;
}

float4x4& operator*=(float4x4& a, const float4x4& b)
{
    a = a * b;
    return a;
}

static float4x4
TranslateTransform(float4 v)
{
    float4x4 ret = Identity4x4;
    for (uint64_t i = 0; i < 4; ++i)
        ((float*)&ret.cols[3])[i] = ((float*)&v)[i];
    
    return ret;
}

static float4x4
ScaleTransform(float4 v)
{
    float4x4 ret = Identity4x4;
    for (uint32_t i = 0; i < 4; ++i)
        ((float*)&ret.cols[i])[i] = ((float*)&v)[i];
    
    return ret;
}

static float4x4
RotateTransform(float radians)
{
    float4x4 ret = Identity4x4;
    
    float c = cosf(radians);
    float s = sinf(radians);

    ((float*)&ret.cols[0])[0] = c;
    ((float*)&ret.cols[0])[1] = s;
    ((float*)&ret.cols[1])[0] = -s;
    ((float*)&ret.cols[1])[1] = c;

    return ret;
}

static float4x4
OrthoProjection(float left, float right, float top, float bottom, float near, float far)
{
    float two = 2.f;
    float rightMinusLeft = (right - left);
    float topMinusBottom = (top - bottom);
    float farMinusNear = (far - near);

    return float4x4 {
        float4{ two / rightMinusLeft, 0.f, 0.f, 0.f},
        float4{0.f, two / topMinusBottom, 0.f, 0.f},
        float4{0.f, 0.f, 1.f / farMinusNear, 0.f},
        float4{-(right + left) / rightMinusLeft, -(top + bottom) / topMinusBottom, -near / farMinusNear, 1.f}
    };
}

static float4x4
PerspectiveProjection(float fovRadians, float aspect, float nearPlane, float farPlane)
{
    float halfFov = fovRadians * .5f;
    float f = 1.f / tanf(halfFov);
    float farMinusNear = farPlane - nearPlane;

    return {
        f / aspect,
        0.f,
        0.f,
        0.f,
        //
        0.f,
        f,
        0.f,
        0.f,
        //
        0.f,
        0.f,
        farPlane / farMinusNear,
        1.f,
        //
        0.f,
        0.f,
        -(farPlane * nearPlane) / farMinusNear,
        0.f,
    };
}


static float4x4
LookAtView(float3 eye, float3 center, float3 up)
{
    float3 f = Norm(center - eye);
    float3 s = Norm(Cross(up, f));
    float3 u = Cross(s, f);

    return {
        s.x,
        u.x,
        f.x,
        0.f,

        s.y,
        u.y,
        f.y,
        0.f,

        s.z,
        u.z,
        f.z,
        0.f,

        -(Dot(s, eye)),
        -(Dot(u, eye)),
        -(Dot(f, eye)),
        1.f,
    };
}

static void
Invert(const float4x4 *m, float4x4 *invOut)
{
    float a[4][8];

    for (int row = 0; row < 4; ++row) {
        for (int col = 0; col < 4; ++col)
            a[row][col] = ((const float *)&m->cols[col])[row];

        for (int col = 0; col < 4; ++col)
            a[row][col + 4] = row == col ? 1.0f : 0.0f;
    }

    for (int col = 0; col < 4; ++col) {
        int pivotRow = col;

        for (int row = col + 1; row < 4; ++row) {
            if (fabsf(a[row][col]) > fabsf(a[pivotRow][col]))
                pivotRow = row;
        }

        ASSERT(fabsf(a[pivotRow][col]) > 1e-8f);

        if (pivotRow != col) {
            for (int i = 0; i < 8; ++i) {
                float temp = a[col][i];
                a[col][i] = a[pivotRow][i];
                a[pivotRow][i] = temp;
            }
        }

        float pivot = a[col][col];

        for (int i = 0; i < 8; ++i)
            a[col][i] /= pivot;

        for (int row = 0; row < 4; ++row) {
            if (row == col)
                continue;

            float factor = a[row][col];

            for (int i = 0; i < 8; ++i)
                a[row][i] -= factor * a[col][i];
        }
    }

    for (int row = 0; row < 4; ++row) {
        for (int col = 0; col < 4; ++col)
            ((float *)&invOut->cols[col])[row] = a[row][col + 4];
    }
}