diff --git a/linalg/Makefile b/linalg/Makefile
new file mode 100644
index 0000000..40b00ed
--- /dev/null
+++ b/linalg/Makefile
@@ -0,0 +1,14 @@
+CC = gcc
+CFLAGS = -Wall -O2
+
+TARGET = main
+SRC = main.c
+
+#LDFLAGS =
+
+$(TARGET): $(SRC)
+	@echo CC $@
+	@$(CC) $(CFLAGS) -o $@ $^ $(LDFLAGS)
+
+clean:
+	rm -f $(TARGET)
diff --git a/linalg/main b/linalg/main
new file mode 100755
index 0000000..7f8b69a
Binary files /dev/null and b/linalg/main differ
diff --git a/linalg/main.c b/linalg/main.c
new file mode 100644
index 0000000..393bf1a
--- /dev/null
+++ b/linalg/main.c
@@ -0,0 +1,238 @@
+#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;
+}