cKit/Foundation/c_Matrix.c

#include <c_Matrix.h>

c_err_t c_Matrix_Init(c_Matrix_t* self, c_size_t rows, c_size_t cols, c_size_t item_size, c_Allocator_t* allocator) {
    if (!self || rows == 0 || cols == 0 || item_size == 0) return C_ERR_PARAM;

    // 防止 rows * cols * item_size 发生整型乘法回绕引发的堆踩踏
    if (rows > (c_size_t)-1 / (cols * item_size)) return C_ERR_OUTOFBOUND;

    self->allocator = (allocator != NULL) ? *allocator : c_DefaultAllocator;
    self->rows = rows;
    self->cols = cols;
    self->item_size = item_size;

    // 一次性开辟完美的连续空间
    c_size_t total_bytes = self->rows * self->cols * self->item_size;
    self->data = c_Allocator_Alloc(&self->allocator, total_bytes);
    if (!self->data) return C_ERR_NOMEM;

    // 默认执行物理抹零清洁
    memset(self->data, 0, total_bytes);
    return C_ERR_OK;
}

// 物理销毁
void c_Matrix_Destroy(c_Matrix_t* self) {
    if (!self) return;
    if (self->data) {
        c_Allocator_Free(&self->allocator, self->data);
        self->data = NULL;
    }
    self->rows = 0;
    self->cols = 0;
}

// 定点安全写入 (值深度复制)
c_err_t c_Matrix_Write(c_Matrix_t* self, c_size_t r, c_size_t c, const void* item) {
    if (!self || !item || r >= self->rows || c >= self->cols) return C_ERR_PARAM;

    // 行优先扁平化寻址计算
    c_size_t index = r * self->cols + c;
    char* target = (char*)self->data + (index * self->item_size);
    memcpy(target, item, self->item_size);

    return C_ERR_OK;
}

// 定点安全读取 (安全拷贝副本)
c_err_t c_Matrix_Read(const c_Matrix_t* self, c_size_t r, c_size_t c, void* out_item) {
    if (!self || !out_item || r >= self->rows || c >= self->cols) return C_ERR_PARAM;

    c_size_t index = r * self->cols + c;
    const char* source = (const char*)self->data + (index * self->item_size);
    memcpy(out_item, source, self->item_size);

    return C_ERR_OK;
}

// 只读原位窥探指针 (由于是固定一维连续空间，只要不销毁，此指针极其平稳安全)
void* c_Matrix_Get(const c_Matrix_t* self, c_size_t r, c_size_t c) {
    if (!self || r >= self->rows || c >= self->cols) return NULL;
    c_size_t index = r * self->cols + c;
    return (char*)self->data + (index * self->item_size);
}

// 批量全局抹值填充
c_err_t c_Matrix_Fill(c_Matrix_t* self, const void* item) {
    if (!self || !item) return C_ERR_PARAM;

    c_size_t total_elements = self->rows * self->cols;
    char* target = (char*)self->data;

    for (c_size_t i = 0; i < total_elements; i++) {
        memcpy(target, item, self->item_size);
        target += self->item_size;
    }
    return C_ERR_OK;
}

// 极致高速的一整行搬迁导出接口 (常用于图像处理的一行扫描线批量搬运)
c_err_t c_Matrix_CopyRow(const c_Matrix_t* self, c_size_t r, void* out_row_buffer) {
    if (!self || !out_row_buffer || r >= self->rows) return C_ERR_PARAM;

    // 因为行优先存储，整行的数据在物理上是 100% 绝对连续的一条线
    // 我们可以直接通过一次 memcpy 瞬时打包带走一整行，效率达到硬件总线传输极限
    const char* row_start = (const char*)self->data + (r * self->cols * self->item_size);
    size_t row_bytes = self->cols * self->item_size;

    memcpy(out_row_buffer, row_start, row_bytes);
    return C_ERR_OK;
}

/* ------------------------------------------------------------------------------------------------------------------ */
/*  */

// 辅助寻址内部宏
#define MAT_ELEMENT(mat, idx) ((char*)(mat)->data + ((idx) * (mat)->item_size))

/**
 * @brief 逐元素（Element-wise）矩阵运算通用驱动函数（内部私有）
 */
static c_err_t c_Matrix_ElementWiseCore(c_Matrix_t* out, const c_Matrix_t* lhs, const c_Matrix_t* rhs,
                                        c_Matrix_OpFn_t op, void* ud) {
    if (!out || !lhs || !rhs || !op) return C_ERR_PARAM;
    // 强御安全校验：矩阵运算要求维度必须完全空间对齐一致
    if (lhs->rows != rhs->rows || lhs->cols != rhs->cols ||
        lhs->rows != out->rows || lhs->cols != out->cols ||
        lhs->item_size != rhs->item_size || lhs->item_size != out->item_size) {
        return C_ERR_OUTOFBOUND;
    }

    c_size_t total_elements = lhs->rows * lhs->cols;

    // 【工业级内核性能极致拉平】：直接抛弃二维坐标计算，用一维指针流极速推进
    for (c_size_t i = 0; i < total_elements; i++) {
        void* out_ptr = MAT_ELEMENT(out, i);
        const void* lhs_ptr = MAT_ELEMENT(lhs, i);
        const void* rhs_ptr = MAT_ELEMENT(rhs, i);
        op(out_ptr, lhs_ptr, rhs_ptr, ud); // 驱动具体类型算子
    }
    return C_ERR_OK;
}

c_err_t c_Matrix_Add(c_Matrix_t* out, const c_Matrix_t* lhs, const c_Matrix_t* rhs, c_Matrix_OpFn_t add_op, void* ud) {
    return c_Matrix_ElementWiseCore(out, lhs, rhs, add_op, ud);
}

c_err_t c_Matrix_Sub(c_Matrix_t* out, const c_Matrix_t* lhs, const c_Matrix_t* rhs, c_Matrix_OpFn_t sub_op, void* ud) {
    return c_Matrix_ElementWiseCore(out, lhs, rhs, sub_op, ud);
}

c_err_t c_Matrix_Div(c_Matrix_t* out, const c_Matrix_t* lhs, const c_Matrix_t* rhs, c_Matrix_OpFn_t div_op, void* ud) {
    return c_Matrix_ElementWiseCore(out, lhs, rhs, div_op, ud);
}

/**
 * @brief 经典线性代数标准矩阵乘法 (Matrix Multiplication: O(N^3))
 * @note 满足拓扑条件：lhs(M x K) * rhs(K x N) = out(M x N)
 */
c_err_t c_Matrix_Mul(c_Matrix_t* out, const c_Matrix_t* lhs, const c_Matrix_t* rhs,
                     c_Matrix_OpFn_t mul_op, c_Matrix_OpFn_t add_op, void* ud) {
    if (!out || !lhs || !rhs || !mul_op || !add_op) return C_ERR_PARAM;
    // 拓扑几何边界校验
    if (lhs->cols != rhs->rows || out->rows != lhs->rows || out->cols != rhs->cols ||
        lhs->item_size != rhs->item_size || lhs->item_size != out->item_size) {
        return C_ERR_OUTOFBOUND;
    }

    c_size_t M = lhs->rows;
    c_size_t K = lhs->cols;
    c_size_t N = rhs->cols;

    // 分配一个栈上的临时缓冲区，用来暂存单次乘法的中间值（避免频繁分配堆内存导致碎片）
    // 由于是泛型，大小通过 item_size 动态抹平
    void* temp_mul_res = c_Allocator_Alloc(&out->allocator, out->item_size);
    if (!temp_mul_res) return C_ERR_NOMEM;

    for (c_size_t i = 0; i < M; i++) {
        for (c_size_t j = 0; j < N; j++) {
            void* out_cell = (char*)out->data + ((i * N + j) * out->item_size);

            // 每次计算新格子前，首先执行干净的物理抹零（清空历史残存值）
            memset(out_cell, 0, out->item_size);

            for (c_size_t k = 0; k < K; k++) {
                const void* lhs_cell = (char*)lhs->data + ((i * K + k) * lhs->item_size);
                const void* rhs_cell = (char*)rhs->data + ((k * N + j) * rhs->item_size);

                // 1. 计算当前权重的乘积：temp_mul_res = lhs_cell * rhs_cell
                mul_op(temp_mul_res, lhs_cell, rhs_cell, ud);

                // 2. 累加到目标格子中：out_cell = out_cell + temp_mul_res
                add_op(out_cell, out_cell, temp_mul_res, ud);
            }
        }
    }

    c_Allocator_Free(&out->allocator, temp_mul_res);
    return C_ERR_OK;
}

/**
 * @brief 标量运算通用驱动核心（内部私有）
 */
C_STATIC_FORCE_INLINE
c_err_t c_Matrix_ScalarCore(c_Matrix_t* out, const c_Matrix_t* in, const void* scalar,
                                   c_Matrix_OpFn_t op, void* ud) {
    if (!out || !in || !scalar || !op) return C_ERR_PARAM;
    if (in->rows != out->rows || in->cols != out->cols || in->item_size != out->item_size) {
        return C_ERR_OUTOFBOUND;
    }

    c_size_t total_elements = in->rows * in->cols;
    for (c_size_t i = 0; i < total_elements; i++) {
        void* out_ptr = MAT_ELEMENT(out, i);
        const void* in_ptr = MAT_ELEMENT(in, i);
        op(out_ptr, in_ptr, scalar, ud); // 驱动单元素与标量运算
    }
    return C_ERR_OK;
}

c_err_t c_Matrix_AddScalar(c_Matrix_t* out, const c_Matrix_t* in, const void* scalar, c_Matrix_OpFn_t add_op, void* ud) {
    return c_Matrix_ScalarCore(out, in, scalar, add_op, ud);
}

c_err_t c_Matrix_SubScalar(c_Matrix_t* out, const c_Matrix_t* in, const void* scalar, c_Matrix_OpFn_t sub_op, void* ud) {
    return c_Matrix_ScalarCore(out, in, scalar, sub_op, ud);
}

c_err_t c_Matrix_MulScalar(c_Matrix_t* out, const c_Matrix_t* in, const void* scalar, c_Matrix_OpFn_t mul_op, void* ud) {
    return c_Matrix_ScalarCore(out, in, scalar, mul_op, ud);
}

c_err_t c_Matrix_DivScalar(c_Matrix_t* out, const c_Matrix_t* in, const void* scalar, c_Matrix_OpFn_t div_op, void* ud) {
    return c_Matrix_ScalarCore(out, in, scalar, div_op, ud);
}

/* ------------------------------------------------------------------------------------------------------------------ */
/*  */

// 矩阵转置实现 (M x N -> N x M)
c_err_t c_Matrix_Transpose(c_Matrix_t* out, const c_Matrix_t* in) {
    if (!out || !in || !out->data || !in->data) return C_ERR_PARAM;

    // 空间几何拓扑强校验
    if (out->rows != in->cols || out->cols != in->rows || out->item_size != in->item_size) {
        return C_ERR_OUTOFBOUND;
    }

    c_size_t r_max = in->rows;
    c_size_t c_max = in->cols;
    c_size_t item_sz = in->item_size;

    // 行优先数据错位映射搬运：in(r, c) -> out(c, r)
    for (c_size_t r = 0; r < r_max; r++) {
        for (c_size_t c = 0; c < c_max; c++) {
            const char* src_cell = (const char*)in->data + ((r * c_max + c) * item_sz);
            char* dst_cell = (char*)out->data + ((c * r_max + r) * item_sz);
            memcpy(dst_cell, src_cell, item_sz); // 泛型值原装迁徙
        }
    }
    return C_ERR_OK;
}

/**
 * @brief 内部私有递归核心：拉普拉斯代数余子式展开
 */
static void c_Matrix_DetInternal(c_Matrix_t* mat, void* out_det,
                                 c_Matrix_OpFn_t add_op, c_Matrix_OpFn_t sub_op,
                                 c_Matrix_OpFn_t mul_op, c_Matrix_NegFn_t neg_op,
                                 void* ud) {
    c_size_t n = mat->rows;
    c_size_t item_sz = mat->item_size;

    // 基本边界 1: 1x1 方阵，行列式值即为其唯一的单体元素本身
    if (n == 1) {
        memcpy(out_det, mat->data, item_sz);
        return;
    }

    // 基本边界 2: 2x2 方阵，计算 ad - bc 快速通路，消除不必要的低层递归
    if (n == 2) {
        char* a = (char*)mat->data + (0 * item_sz);
        char* b = (char*)mat->data + (1 * item_sz);
        char* c = (char*)mat->data + (2 * item_sz);
        char* d = (char*)mat->data + (3 * item_sz);

        // 分配栈上局部存储，避开堆碎片
        void* ad = c_Allocator_Alloc(&mat->allocator, item_sz);
        void* bc = c_Allocator_Alloc(&mat->allocator, item_sz);

        mul_op(ad, a, d, ud); // a * d
        mul_op(bc, b, c, ud); // b * c
        sub_op(out_det, ad, bc, ud); // ad - bc

        c_Allocator_Free(&mat->allocator, ad);
        c_Allocator_Free(&mat->allocator, bc);
        return;
    }

    // 递归分支: n > 2，固定沿第 0 行进行展开
    memset(out_det, 0, item_sz); // 抹零累加器

    // 原地创建子方阵控制头 (大小为 n-1)
    c_Matrix_t sub_mat;
    sub_mat.allocator = mat->allocator;
    sub_mat.rows = n - 1;
    sub_mat.cols = n - 1;
    sub_mat.item_size = item_sz;

    c_size_t sub_bytes = sub_mat.rows * sub_mat.cols * item_sz;
    sub_mat.data = c_Allocator_Alloc(&mat->allocator, sub_bytes);
    if (!sub_mat.data) return;

    void* sub_det = c_Allocator_Alloc(&mat->allocator, item_sz);
    void* term = c_Allocator_Alloc(&mat->allocator, item_sz);
    void* neg_term = c_Allocator_Alloc(&mat->allocator, item_sz);

    // 遍历第 0 行的每一列 c
    for (c_size_t c = 0; c < n; c++) {
        // 构建割裂后的代数余子式子方阵数据块
        c_size_t sub_r = 0;
        for (c_size_t r = 1; r < n; r++) { // 跳过第 0 行
            c_size_t sub_c = 0;
            for (c_size_t j = 0; j < n; j++) {
                if (j == c) continue; // 跳过当前列

                const char* src = (const char*)mat->data + ((r * n + j) * item_sz);
                char* dst = (char*)sub_mat.data + ((sub_r * (n - 1) + sub_c) * item_sz);
                memcpy(dst, src, item_sz);
                sub_c++;
            }
            sub_r++;
        }

        // 递归求解子方阵的行列式
        c_Matrix_DetInternal(&sub_mat, sub_det, add_op, sub_op, mul_op, neg_op, ud);

        // 计算当前项的系数乘积 term = mat(0, c) * sub_det
        const void* current_element = (const char*)mat->data + (c * item_sz);
        mul_op(term, current_element, sub_det, ud);

        // 根据棋盘格正负号规则 ((-1)^(r+c))：当前第 0 行第 c 列，当 c 为奇数时取反
        if (c % 2 == 1) {
            neg_op(neg_term, term, ud); // 取负号
            add_op(out_det, out_det, neg_term, ud); // 累加负项
        } else {
            add_op(out_det, out_det, term, ud); // 累加正项
        }
    }

    // 异常安全性层级强力物理释放，绝不产生多级残留
    c_Allocator_Free(&mat->allocator, sub_mat.data);
    c_Allocator_Free(&mat->allocator, sub_det);
    c_Allocator_Free(&mat->allocator, term);
    c_Allocator_Free(&mat->allocator, neg_term);
}

// 行列式主驱动包装接口
c_err_t c_Matrix_Determinant(c_Matrix_t* self, void* out_det,
                             c_Matrix_OpFn_t add_op, c_Matrix_OpFn_t sub_op,
                             c_Matrix_OpFn_t mul_op, c_Matrix_NegFn_t neg_op,
                             void* ud) {
    if (!self || !out_det || !add_op || !sub_op || !mul_op || !neg_op) return C_ERR_PARAM;

    // 强行约束约束：只有正方矩阵具备行列式
    if (self->rows != self->cols || self->rows == 0) {
        return C_ERR_OUTOFBOUND;
    }

    c_Matrix_DetInternal(self, out_det, add_op, sub_op, mul_op, neg_op, ud);
    return C_ERR_OK;
}