#include <assert.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>

/**
 * @brief A 2×2 matrix stored in row-major order.
 *
 * The elements are laid out as:
 * [ arr[0] arr[1] ]
 * [ arr[2] arr[3] ]
 */
typedef struct Mat2 {
    float arr[4];
} Mat2;

/**
 * @brief A 3×3 matrix stored in row-major order.
 *
 * The elements are laid out as:
 * [ arr[0] arr[1] arr[2] ]
 * [ arr[3] arr[4] arr[5] ]
 * [ arr[6] arr[7] arr[8] ]
 */
typedef struct Mat3 {
    float arr[9];
} Mat3;

/**
 * @brief A 2D vector with x and y components.
 */
typedef struct Vec2 {
    float x, y;
} Vec2;

/**
 * @brief A 3D vector with x, y, and z components.
 */
typedef struct Vec3 {
    float x, y, z;
} Vec3;

/**
 * @brief Computes the dot product of two 3D vectors.
 *
 * @param a Pointer to the first vector.
 * @param b Pointer to the second vector.
 * @return The dot product (a • b).
 */
inline float vec3_dot(const Vec3 *a, const Vec3 *b);

/**
 * @brief Computes the cross product of two 3D vectors.
 *
 * @param a Pointer to the first vector.
 * @param b Pointer to the second vector.
 * @return The cross product vector (a × b).
 */
Vec3 vec3_cross(const Vec3 *a, const Vec3 *b);

/**
 * @brief Computes the determinant of a 2×2 matrix.
 *
 * @param m Pointer to the matrix.
 * @return The determinant of the matrix.
 */
float mat2_det(const Mat2 *m);

/**
 * @brief Computes the determinant of a 3×3 matrix (row-major order).
 *
 * @param m Pointer to the matrix.
 * @return The determinant of the matrix.
 */
float mat3_det2(const Mat3 *m);

/**
 * @brief Multiplies two 2×2 matrices (row-major order).
 *
 * @param m1 Pointer to the first matrix.
 * @param m2 Pointer to the second matrix.
 * @return The resulting matrix product (m1 × m2).
 */
Mat2 mat2_mul(const Mat2 *m1, const Mat2 *m2);

/**
 * @brief Multiplies two 3×3 matrices (row-major order).
 *
 * @param m1 Pointer to the first matrix.
 * @param m2 Pointer to the second matrix.
 * @return The resulting matrix product (m1 × m2).
 */
Mat3 mat3_mul(const Mat3 *m1, const Mat3 *m2);

/**
 * @brief Checks if two 2×2 matrices are approximately equal.
 *
 * @param a Pointer to the first matrix.
 * @param b Pointer to the second matrix.
 * @param epsilon Tolerance for comparison.
 * @return true if all elements are approximately equal within epsilon.
 */
bool mat2_approx_eq(const Mat2 *a, const Mat2 *b, float epsilon);

/**
 * @brief Checks if two 3×3 matrices are approximately equal.
 *
 * @param a Pointer to the first matrix.
 * @param b Pointer to the second matrix.
 * @param epsilon Tolerance for comparison.
 * @return true if all elements are approximately equal within epsilon.
 */
bool mat3_approx_eq(const Mat3 *a, const Mat3 *b, float epsilon);

#define MAT2_AT(m, row, col) ((m)->arr[(col) * 2 + (row)])
#define MAT3_AT(m, row, col) ((m)->arr[(col) * 3 + (row)])

/* Header end... */

float mat3_det(const Mat3 *m) {
    float m00 = MAT3_AT(m, 0, 0);
    float m01 = MAT3_AT(m, 0, 1);
    float m02 = MAT3_AT(m, 0, 2);

    float m10 = MAT3_AT(m, 1, 0);
    float m11 = MAT3_AT(m, 1, 1);
    float m12 = MAT3_AT(m, 1, 2);

    float m20 = MAT3_AT(m, 2, 0);
    float m21 = MAT3_AT(m, 2, 1);
    float m22 = MAT3_AT(m, 2, 2);

    return m00 * (m11 * m22 - m12 * m21) - m01 * (m10 * m22 - m12 * m20) +
           m02 * (m10 * m21 - m11 * m20);
}

inline float vec3_dot(const struct Vec3 *a, const struct Vec3 *b) {
    return a->x * b->x + a->y * b->y + a->z * b->z;
}

struct Vec3 vec3_cross(const struct Vec3 *a, const struct Vec3 *b) {
    struct Vec3 res = {.x = a->y * b->z - a->z * b->y,
                       .y = a->x * b->z - a->z * b->x,
                       .z = a->x * b->y - a->y * b->x};

    return res;
}

Mat2 mat2_mul(const Mat2 *m1, const Mat2 *m2) {
    Mat2 m3 = {.arr = {
                   MAT2_AT(m1, 0, 0) * MAT2_AT(m2, 0, 0) +
                       MAT2_AT(m1, 0, 1) * MAT2_AT(m2, 1, 0),
                   MAT2_AT(m1, 1, 0) * MAT2_AT(m2, 0, 0) +
                       MAT2_AT(m1, 1, 1) * MAT2_AT(m2, 1, 0),
                   MAT2_AT(m1, 0, 0) * MAT2_AT(m2, 0, 1) +
                       MAT2_AT(m1, 0, 1) * MAT2_AT(m2, 1, 1),
                   MAT2_AT(m1, 1, 0) * MAT2_AT(m2, 0, 1) +
                       MAT2_AT(m1, 1, 1) * MAT2_AT(m2, 1, 1),
               }};

    return m3;
}

Mat3 mat3_mul(const Mat3 *m1, const Mat3 *m2) {
    Mat3 m3;

    for (int col = 0; col < 3; ++col) {
        for (int row = 0; row < 3; ++row) {
            float sum = 0.0f;
            for (int k = 0; k < 3; ++k) {
                sum += MAT3_AT(m1, row, k) * MAT3_AT(m2, k, col);
            }
            m3.arr[col * 3 + row] = sum;
        }
    }

    return m3;
}

void mat3_print(const Mat3 *m) {
    for (int row = 0; row < 3; ++row) {
        printf("| ");
        for (int col = 0; col < 3; ++col) {
            printf("%8.3f ", MAT3_AT(m, row, col));
        }
        printf("|\n");
    }
}

bool mat2_approx_eq(const Mat2 *a, const Mat2 *b, float epsilon) {
    for (int i = 0; i < 4; ++i) {
        if (fabsf(a->arr[i] - b->arr[i]) > epsilon)
            return false;
    }
    return true;
}

bool mat3_approx_eq(const Mat3 *a, const Mat3 *b, float epsilon) {
    for (int i = 0; i < 9; ++i) {
        if (fabsf(a->arr[i] - b->arr[i]) > epsilon)
            return false;
    }
    return true;
}

/* Implem end */

int main(void) {
    {
        Mat3 m = {{1, 0, 0, 0, 1, 0, 0, 0, 1}};

        Mat3 m3 = mat3_mul(&m, &m);
        assert(mat3_approx_eq(&m3, &m, 0.01));
    }
    {
        Mat3 m = {{1, 0, 0, 0, 1, 0, 0, 0, 1}};
        float d = mat3_det(&m);
        printf("Determinant: %f\n", d);

        MAT3_AT(&m, 0, 0) = 2;
        d = mat3_det(&m);
        printf("Determinant: %f\n", d);
    }
    {
        struct Vec3 a = {10, 10, 10};
        struct Vec3 b = {5, 5, 5};
        struct Vec3 c = vec3_cross(&a, &b);
        printf("{ Vec3: %f, %f, %f }\n", c.x, c.y, c.z);
    }
    {
        struct Vec3 a = {0, 1, 0};
        struct Vec3 b = {0, 0, 1};
        struct Vec3 c = vec3_cross(&a, &b);
        printf("{ Vec3: %f, %f, %f }\n", c.x, c.y, c.z);
    }

    return 0;
}