增加了指定nvcc路径以及主机编译器的指令

.gitignore

@@ -1 +1,2 @@
-build/
+build/
+.vscode/

CMakeLists.txt

@@ -1,27 +1,33 @@
 cmake_minimum_required(VERSION 3.24)  # 设置cmake最低版本
 
+# 设置 nvcc 的路径
+set(CMAKE_CUDA_COMPILER /usr/local/cuda/bin/nvcc)  # 根据实际安装路径设置
+
+# 设置 nvcc 使用的主机编译器
+set(CMAKE_CUDA_HOST_COMPILER /usr/bin/g++-11)
+
 # 设置 C++ 标准
-set(CMAKE_CXX_STANDARD 20)
-set(CMAKE_CXX_STANDARD_REQUIRED True)
+# set(CMAKE_CXX_STANDARD 20)
+# set(CMAKE_CXX_STANDARD_REQUIRED True)
 # 设置 CUDA C++ 标准
 set(CMAKE_CUDA_STANDARD 20)
 set(CMAKE_CUDA_STANDARD_REQUIRED True)
 set(CMAKE_CUDA_ARCHITECTURES native)  # 设置CUDA架构
 
 # 设置编译选项
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
-set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3")
+# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
+set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 -maxrregcount=128")
 
-project(cuElim_GF28 LANGUAGES CXX CUDA)  # 设置项目名和语言
+project(cuElim LANGUAGES CXX CUDA)  # 设置项目名和语言
 
 include_directories(${PROJECT_SOURCE_DIR}/include)  # 添加头文件目录
 
 find_package(OpenMP REQUIRED) # 查找 OpenMP 包
 
 add_executable(cuelim ./src/main.cu)  # 添加可执行文件
-target_link_libraries(cuelim OpenMP::OpenMP_CXX)  # 链接 OpenMP 库
+target_link_libraries(cuelim OpenMP::OpenMP_CXX m4ri)  # 链接 OpenMP 库
 
 enable_testing()  # 启动ctest测试
 add_subdirectory(test)  # 添加测试目录
 
-# add_subdirectory(benchmark)  # 添加性能测试目录
+add_subdirectory(benchmark)  # 添加性能测试目录

README.md

@@ -0,0 +1,52 @@
+# cuElim
+
+使用CUDA统一内存实现GF2域，GF256域和素域上的矩阵乘法和高斯消元，无需考虑显存大小。
+
+## 使用说明
+
+### 加入到现有项目
+
+1. 将`include`文件夹中的所有头文件添加到现有项目中。
+
+2. 在需要使用gpu函数的地方引用`cuelim.cuh`或`cuelim_m4ri.cuh`（GF2域上高斯消元接口）或`cuelim_m4rie.cuh`（GF256域上高斯消元接口）。
+
+### 使用当前项目
+
+1. 安装依赖
+
+   1. `C++ CUDA CMake ...`
+   2. [`GoogleTest`](https://github.com/google/googletest)
+   3. [`GoogleBenchmark`](https://github.com/google/benchmark)
+
+2. 构建项目
+
+   ```sh
+   mkdir build && cd build
+   cmake ..
+   make -j  # 同时编译多个目标
+   ctest  # 或make test 执行所有测试
+   ```
+
+3. 运行可执行文件
+
+   ```sh
+   ./cuelim  # 执行主程序
+   ./test/target  # 执行特定测试
+   ./benchmark/target  # 执行特定性能测试
+   ```
+
+## 功能简介
+
+- `ElimResult`：存储高斯消元的结果，包含秩、主元行（进行行交换前的位置）、主元列。
+
+- `gf2::MatGF2`：储存GF2矩阵，每个uint64_t储存64个元素，单个uint64_t中元素从低位到高位排列。
+- `gf2::ElimResult gf2::MatGF2::gpu_elim()`：高斯消元。
+
+- `gf256::MatGF256`：存储GF256矩阵，每个uint64_t储存8个元素，单个uint64_t中元素从低位到高位排列。
+- `gf256::ElimResult gf256::MatGF256::gpu_elim(const gf256::GF256 &gf)`：进行高斯消元
+
+- `gfp::MatGFP`：储存GF65521矩阵，每个uint32_t储存1个元素
+- `gfp::ElimResult gfp::MatGFP::gpu_elim()`：进行高斯消元
+
+- `size_t gpu_mzd_elim(mzd_t *A)`：m4ri接口
+- `size_t gpu_mzed_elim(mzed_t *A)`：m4rie接口

benchmark/CMakeLists.txt

@@ -0,0 +1,14 @@
+find_package(benchmark REQUIRED)
+
+include_directories(${PROJECT_SOURCE_DIR}/test)
+
+set(BENCH_SRC_FILES
+    "bench_gf256_mul.cu"
+    "bench_gfp_mul.cu"
+)
+
+foreach(SRC ${BENCH_SRC_FILES})
+    get_filename_component(SRC_NAME ${SRC} NAME_WE)
+    add_executable(${SRC_NAME} ${SRC})
+    target_link_libraries(${SRC_NAME} benchmark::benchmark benchmark::benchmark_main)
+endforeach()

benchmark/bench_gf256_mul.cu

@@ -0,0 +1,21 @@
+#include <benchmark/benchmark.h>
+#include "cuelim.cuh"
+
+using namespace gf256;
+
+template <MatGF256 (*GpuFunc)(const MatGF256 &, const MatGF256 &, const GF256 &)>
+void bench_gf256_mul(benchmark::State &state)
+{
+    GF256 ff(0b100011101);
+    uint_fast32_t seed = 41921095;
+    size_t x = state.range(0), y = state.range(1), z = state.range(2);
+    MatGF256 A(x, y), B(y, z);
+    A.randomize(seed);
+    B.randomize(seed);
+    for (auto _ : state)
+    {
+        MatGF256 C(GpuFunc(A, B, ff));
+    }
+}
+
+BENCHMARK_TEMPLATE(bench_gf256_mul, gpu_mul)->Args({10000, 10000, 10000});

benchmark/bench_gfp_mul.cu

@@ -0,0 +1,19 @@
+#include <benchmark/benchmark.h>
+#include "test_header.cuh"
+
+using namespace gfp;
+
+static void bench_gfp(benchmark::State &state)
+{
+    uint_fast32_t seed = 41921095;
+    size_t x = state.range(0), y = state.range(1), z = state.range(2);
+    MatGFP A(x, y), B(y, z);
+    A.randomize(seed);
+    B.randomize(seed);
+    for (auto _ : state)
+    {
+        MatGFP C = A * B;
+    }
+}
+
+BENCHMARK(bench_gfp)->Args({10000, 10000, 10000});

include/cpp_progress.hpp

@@ -0,0 +1,72 @@
+#ifndef CPP_PROGRESS_HPP
+#define CPP_PROGRESS_HPP
+
+#include <iostream>
+#include <cassert>
+#include <cstdint>
+
+namespace progress
+{
+    class ProgressBar
+    {
+    public:
+        ProgressBar(const std::string &desc, const int64_t total_ticks, const int64_t bar_width = 50, const int64_t ticks_per_display = 1) : desc{desc}, total_ticks{total_ticks}, bar_width{bar_width}, ticks_per_display{ticks_per_display}
+        {
+            assert(total_ticks > 0 && bar_width > 0 && ticks_per_display > 0);
+        }
+
+        void tick_display()
+        {
+#ifdef SHOW_PROGRESS_BAR
+            if (++ticks == total_ticks)
+            {
+                done();
+                return;
+            }
+            double progress = static_cast<double>(ticks) / total_ticks;
+            int64_t pos = static_cast<int64_t>(bar_width * progress);
+            display(pos);
+#endif
+        }
+
+    private:
+        std::string get_bar(const int64_t pos)
+        {
+            if (bar != "" && pos == now_pos)
+                return bar;
+            bar.clear();
+            for (int i = 0; i < bar_width; ++i)
+            {
+                if (i < pos)
+                    bar += '=';
+                else if (i == pos)
+                    bar += ">";
+                else
+                    bar += ' ';
+            }
+            now_pos = pos;
+            return bar;
+        }
+
+        void display(int64_t pos)
+        {
+            std::cout << "\33[2K\r" << "[" << get_bar(pos) << "]" << desc << std::flush;
+        }
+
+        void done()
+        {
+            display(bar_width);
+            std::cout << std::endl;
+        }
+
+        int64_t ticks = 0;
+        int64_t now_pos = -1;
+        std::string bar = "";
+        std::string desc = "";
+        const int64_t total_ticks;
+        const int64_t bar_width;
+        int64_t ticks_per_display;
+    };
+}
+
+#endif

include/cuelim.cuh

@@ -1,6 +1,13 @@
 #ifndef CUELIM_CUH
 #define CUELIM_CUH
 
-#include "multiplication.cuh"
+#include "gf2/gf2_mul.cuh"
+#include "gf2/gf2_elim.cuh"
+
+#include "gf256/gf256_mul.cuh"
+#include "gf256/gf256_elim.cuh"
+
+#include "gfp/gfp_mul.cuh"
+#include "gfp/gfp_elim.cuh"
 
 #endif

include/cuelim_m4ri.cuh

@@ -0,0 +1,44 @@
+#ifndef CUELIM_M4RI_CUH
+#define CUELIM_M4RI_CUH
+
+#include "gf2/gf2_mul.cuh"
+#include "gf2/gf2_elim.cuh"
+#include <m4ri/m4ri.h>
+
+namespace gf2
+{
+    void mzdread(mzd_t *A, MatGF2 &mat)
+    {
+        assert(A->nrows == mat.nrows && A->ncols == mat.ncols);
+        for (size_t r = 0; r < A->nrows; r++)
+        {
+            for (size_t cn = 0; cn < A->width; cn++)
+            {
+                *mat.at_base(r, cn) = A->rows[r][cn];
+            }
+        }
+    }
+
+    void mzdwrite(MatGF2 &mat, mzd_t *A)
+    {
+        assert(A->nrows == mat.nrows && A->ncols == mat.ncols);
+        for (size_t r = 0; r < mat.nrows; r++)
+        {
+            for (size_t cn = 0; cn < mat.width; cn++)
+            {
+                A->rows[r][cn] = *mat.at_base(r, cn);
+            }
+        }
+    }
+}
+
+size_t gpu_mzd_elim(mzd_t *A)
+{
+    gf2::MatGF2 mat(A->nrows, A->ncols);
+    gf2::mzdread(A, mat);
+    gf2::ElimResult res = mat.gpu_elim();
+    gf2::mzdwrite(mat, A);
+    return res.rank;
+}
+
+#endif

include/cuelim_m4rie.cuh

@@ -0,0 +1,45 @@
+#ifndef INTERFACE_CUH
+#define INTERFACE_CUH
+
+#include "gf256/gf256_mul.cuh"
+#include "gf256/gf256_elim.cuh"
+#include <m4rie/m4rie.h>
+
+namespace gf256
+{
+    void mzedread(mzed_t *A, MatGF256 &mat)
+    {
+        assert(A->nrows == mat.nrows && A->ncols == mat.ncols);
+        for (size_t r = 0; r < A->nrows; r++)
+        {
+            for (size_t cn = 0; cn < A->x->width; cn++)
+            {
+                *mat.at_base(r, cn) = A->x->rows[r][cn];
+            }
+        }
+    }
+
+    void mzedwrite(MatGF256 &mat, mzed_t *A)
+    {
+        assert(A->nrows == mat.nrows && A->ncols == mat.ncols);
+        for (size_t r = 0; r < mat.nrows; r++)
+        {
+            for (size_t cn = 0; cn < mat.width; cn++)
+            {
+                A->x->rows[r][cn] = *mat.at_base(r, cn);
+            }
+        }
+    }
+}
+
+size_t gpu_mzed_elim(mzed_t *A)
+{
+    gf256::MatGF256 mat(A->nrows, A->ncols);
+    gf256::mzedread(A, mat);
+    gf256::GF256 gf256(A->finite_field->minpoly);
+    gf256::ElimResult res = mat.gpu_elim(gf256);
+    gf256::mzedwrite(mat, A);
+    return res.rank;
+}
+
+#endif

include/gf2/gf2_elim.cuh

@@ -0,0 +1,207 @@
+#ifndef GF2_ELIM_CUH
+#define GF2_ELIM_CUH
+
+#include "gf2_mat.cuh"
+
+namespace gf2
+{
+    size_t cpu_elim_base(base_t *base_col, base_t base_col_len, size_t st_r, size_t w, vector<size_t> &p_col, vector<size_t> &p_row)
+    {
+        size_t rank = 0;
+        size_t pivot[gf2_num];
+        size_t next[gf2_num];
+        for (size_t pivot_col = 0; pivot_col < gf2_num; pivot_col++)
+        {
+            for (size_t r = rank; r < base_col_len; r++)
+            {
+                for (size_t i = 0; i < rank; i++)
+                {
+                    if (next[i] == r)
+                    {
+                        if (get(base_col[r], pivot[i]) != 0)
+                        {
+                            base_col[r] ^= concat(base_zero, pivot[i] + 1, base_col[i]);
+                        }
+                        next[i]++;
+                    }
+                }
+
+                if (get(base_col[r], pivot_col) != 0)
+                {
+                    p_col.push_back(w * gf2_num + pivot_col);
+                    p_row.push_back(st_r + r);
+                    if (r != rank)
+                    {
+                        base_t temp = base_col[rank];
+                        base_col[rank] = base_col[r];
+                        base_col[r] = temp;
+                    }
+                    pivot[rank] = pivot_col;
+                    next[rank] = rank + 1;
+                    rank++;
+                    break;
+                }
+            }
+        }
+        return rank;
+    }
+
+    __managed__ uint32_t m_pivot[gf2_num];
+
+    __global__ void gpu_mksrc_kernel(base_t *src, size_t s_rowstride, base_t *base_col, size_t rank, uint32_t m_pivot[gf2_num], size_t width)
+    {
+        size_t w = blockIdx.x * blockDim.x + threadIdx.x;
+        if (w >= width)
+        {
+            return;
+        }
+        base_t temp[gf2_num];
+        for (size_t r = 0; r < rank; r++)
+        {
+            temp[r] = *at_base(src, s_rowstride, r, w);
+        }
+        for (size_t r = 0; r < rank; r++)
+        {
+            for (size_t i = 0; i < r; i++)
+            {
+                if (get(base_col[r], m_pivot[i]) != 0)
+                {
+                    temp[r] ^= temp[i];
+                }
+            }
+        }
+        for (size_t rr = 1; rr < rank; rr++)
+        {
+            size_t r = rank - 1 - rr;
+            for (size_t i = r + 1; i < rank; i++)
+            {
+                if (get(base_col[r], m_pivot[i]) != 0)
+                {
+                    temp[r] ^= temp[i];
+                }
+            }
+        }
+        for (size_t r = 0; r < rank; r++)
+        {
+            *at_base(src, s_rowstride, r, w) = temp[r];
+        }
+    }
+
+    __global__ void gpu_elim_mktb_kernel(base_t *tb, size_t tb_rowstride, base_t *b, size_t b_rowstride, size_t tb_width)
+    {
+        size_t w = blockIdx.x * blockDim.x + threadIdx.x;
+        size_t r = blockIdx.y * blockDim.y + threadIdx.y;
+
+        if (w >= tb_width)
+        {
+            return;
+        }
+
+        base_t val = base_zero;
+        base_t idx = r & gf2_table_mask;
+        base_t st_row = (r >> gf2_table_len) * gf2_table_len;
+
+        for (size_t i = 0; i < gf2_table_len; i++)
+        {
+            if (get(idx, i) != 0)
+            {
+                val ^= *at_base(b, b_rowstride, st_row + i, w);
+            }
+        }
+        *at_base(tb, tb_rowstride, r, w) = val;
+    }
+
+    __global__ void gpu_elim_kernel(base_t *idx, base_t *tb, size_t tb_rowstride, base_t *data, size_t rowstride, size_t rank, uint32_t m_pivot[gf2_num], size_t st_skip, size_t width, size_t nrows)
+    {
+        size_t w = blockIdx.x * blockDim.x + threadIdx.x;
+        size_t r = blockIdx.y * blockDim.y + threadIdx.y;
+
+        if (w >= width || r >= nrows || (r >= st_skip && r < st_skip + rank))
+        {
+            return;
+        }
+
+        base_t val = idx[r];
+        base_t temp = base_zero;
+        for (size_t i = 0; i < gf2_table_num; i++)
+        {
+            base_t loc = base_zero;
+            for (size_t j = 0; (j < gf2_table_len) && (i * gf2_table_len + j < rank); j++)
+            {
+                loc |= get(val, m_pivot[i * gf2_table_len + j]) << j;
+            }
+            temp ^= *at_base(tb, tb_rowstride, i * (1 << gf2_table_len) + loc, w);
+        }
+        *at_base(data, rowstride, r, w) ^= temp;
+    }
+
+    // __managed__ base_t spL[gf2_num];
+
+    __host__ ElimResult MatGF2::gpu_elim()
+    {
+        MatGF2 tb(gf2_table_num * (1 << gf2_table_len), ncols);
+
+        base_t *base_col;
+        cudaMallocManaged(&base_col, nrows * sizeof(base_t));
+        base_t *idx;
+        cudaMallocManaged(&idx, nrows * sizeof(base_t));
+
+        size_t rank = 0;
+        vector<size_t> p_col, p_row;
+
+        progress::ProgressBar pb("GPU ELIMINATE", width);
+        for (size_t w = 0; w < width; w++, pb.tick_display())
+        {
+            CUDA_CHECK(cudaMemcpy2D(base_col + rank, sizeof(base_t), at_base(rank, w), rowstride * sizeof(base_t), sizeof(base_t), nrows - rank, cudaMemcpyDefault));
+
+            size_t src_rank = cpu_elim_base(base_col + rank, nrows - rank, rank, w, p_col, p_row);
+
+            if (src_rank == 0)
+            {
+                continue;
+            }
+
+            for (size_t i = 0; i < src_rank; i++)
+            {
+                cpu_swap_row(rank + i, p_row[rank + i]);
+            }
+
+            for (size_t r = 0; r < src_rank; r++)
+            {
+                size_t loc = (p_col[rank + r] - w * gf2_num);
+                m_pivot[r] = loc;
+            }
+
+            dim3 block_src(THREAD_X);
+            dim3 grid_src((width - w - 1) / block_src.x + 1);
+            gpu_mksrc_kernel<<<grid_src, block_src>>>(at_base(rank, w), rowstride, base_col + rank, src_rank, m_pivot, width);
+            cudaDeviceSynchronize();
+
+            size_t tb_nrows = (src_rank / gf2_table_len) * (1 << gf2_table_len) + (src_rank % gf2_table_len == 0 ? 0 : 1 << (src_rank % gf2_table_len));
+
+            dim3 block_tb(THREAD_X, THREAD_Y);
+            dim3 grid_tb((width - w - 1) / block_tb.x + 1, (tb_nrows - 1) / block_tb.y + 1);
+            gpu_elim_mktb_kernel<<<grid_tb, block_tb>>>(tb.data, tb.rowstride, at_base(rank, w), rowstride, tb.width);
+            cudaDeviceSynchronize();
+
+            CUDA_CHECK(cudaMemcpy2D(idx, sizeof(base_t), at_base(0, w), rowstride * sizeof(base_t), sizeof(base_t), nrows, cudaMemcpyDefault));
+
+            dim3 block(THREAD_X, THREAD_Y);
+            dim3 grid((width - w - 1) / block.x + 1, (nrows - 1) / block.y + 1);
+            gpu_elim_kernel<<<grid, block>>>(idx, tb.data, tb.rowstride, at_base(0, w), rowstride, src_rank, m_pivot, rank, width - w, nrows);
+            cudaDeviceSynchronize();
+
+            rank += src_rank;
+
+            if (rank == nrows)
+            {
+                break;
+            }
+        }
+        cudaFree(base_col);
+        cudaFree(idx);
+        return {rank, p_col, p_row};
+    }
+}
+
+#endif

include/gf2/gf2_header.cuh

@@ -0,0 +1,51 @@
+#ifndef GF2_HEADER_CUH
+#define GF2_HEADER_CUH
+
+#include "../header.cuh"
+
+namespace gf2
+{
+    static const size_t gf2_num = base_len;
+    static const size_t gf2_len = 1;
+
+    static const size_t gf2_table_num = 8;
+    static const size_t gf2_table_len = 8;
+    static const size_t gf2_table_mask = 0xFF;
+
+    __host__ inline void del(base_t &dst, size_t idx)
+    {
+        dst &= ~(base_one << idx);
+    }
+    __host__ inline void set(base_t &dst, size_t idx)
+    {
+        dst |= base_one << idx;
+    }
+    __host__ inline base_t concat(base_t dst_l, size_t idx_l, base_t dst_r)
+    {
+        if (idx_l == 0)
+        {
+            return dst_r;
+        }
+        if (idx_l == gf2_num)
+        {
+            return dst_l;
+        }
+        return (dst_l & (base_fullmask >> (base_len - idx_l))) | (dst_r & (base_fullmask << idx_l));
+    }
+    __host__ __device__ inline base_t get(base_t src, size_t idx)
+    {
+        return (src >> idx) & base_one;
+    }
+
+    __host__ inline base_t rev(base_t n)
+    {
+        n = (n & (base_t)0xAA'AA'AA'AA'AA'AA'AA'AA) >> 1 | (n & (base_t)0x55'55'55'55'55'55'55'55) << 1;
+        n = (n & (base_t)0xCC'CC'CC'CC'CC'CC'CC'CC) >> 2 | (n & (base_t)0x33'33'33'33'33'33'33'33) << 2;
+        n = (n & (base_t)0xF0'F0'F0'F0'F0'F0'F0'F0) >> 4 | (n & (base_t)0x0F'0F'0F'0F'0F'0F'0F'0F) << 4;
+        n = (n & (base_t)0xFF'00'FF'00'FF'00'FF'00) >> 8 | (n & (base_t)0x00'FF'00'FF'00'FF'00'FF) << 8;
+        n = (n & (base_t)0xFF'FF'00'00'FF'FF'00'00) >> 16 | (n & (base_t)0x00'00'FF'FF'00'00'FF'FF) << 16;
+        return n >> 32 | n << 32;
+    }
+
+}
+#endif

include/gf2/gf2_mat.cuh

@@ -0,0 +1,249 @@
+#ifndef GF2_MAT_CUH
+#define GF2_MAT_CUH
+
+#include "gf2_header.cuh"
+
+#include <random>
+#include <vector>
+#include <bitset>
+// #include <algorithm>
+
+namespace gf2
+{
+    struct ElimResult
+    {
+        size_t rank;
+        vector<size_t> pivot;
+        vector<size_t> swap_row;
+    };
+
+    class MatGF2
+    {
+    public:
+        enum MatType
+        {
+            root,
+            window,
+            moved,
+        };
+        // 只能构造root矩阵
+        MatGF2(size_t nrows, size_t ncols) : nrows(nrows), ncols(ncols), type(root)
+        {
+            width = (ncols - 1) / gf2_num + 1;
+            rowstride = ((width - 1) / 4 + 1) * 4; // 以32字节（4*64bit）为单位对齐
+            CUDA_CHECK(cudaMallocManaged((void **)&data, nrows * rowstride * sizeof(base_t)));
+            CUDA_CHECK(cudaMemset(data, 0, nrows * rowstride * sizeof(base_t)));
+        }
+        // 只能以base_t为单位建立window矩阵
+        MatGF2(const MatGF2 &src, size_t begin_ri, size_t begin_wi, size_t end_rj, size_t end_wj) : nrows(end_rj - begin_ri), ncols((end_wj == src.width ? src.ncols : end_wj * gf2_num) - begin_wi * gf2_num), width(end_wj - begin_wi), rowstride(src.rowstride), type(window), data(src.at_base(begin_ri, begin_wi))
+        {
+            assert(begin_ri < end_rj && end_rj <= src.nrows && begin_wi < end_wj && end_wj <= src.width);
+        }
+        // 只能拷贝构造root矩阵
+        MatGF2(const MatGF2 &m) : MatGF2(m.nrows, m.ncols)
+        {
+            CUDA_CHECK(cudaMemcpy2D(data, rowstride * sizeof(base_t), m.data, m.rowstride * sizeof(base_t), m.width * sizeof(base_t), nrows, cudaMemcpyDefault));
+        }
+        MatGF2(MatGF2 &&m) noexcept : nrows(m.nrows), ncols(m.ncols), width(m.width), rowstride(m.rowstride), type(m.type), data(m.data)
+        {
+            m.type = moved;
+            m.data = nullptr;
+        }
+        MatGF2 &operator=(const MatGF2 &m)
+        {
+            if (this == &m)
+            {
+                return *this;
+            }
+            assert(nrows == m.nrows && ncols == m.ncols);
+            CUDA_CHECK(cudaMemcpy2D(data, rowstride * sizeof(base_t), m.data, m.rowstride * sizeof(base_t), m.width * sizeof(base_t), nrows, cudaMemcpyDefault));
+            return *this;
+        }
+        MatGF2 &operator=(MatGF2 &&m) noexcept
+        {
+            if (this == &m)
+            {
+                return *this;
+            }
+            if (type == root)
+            {
+                CUDA_CHECK(cudaFree(data));
+            }
+            nrows = m.nrows;
+            ncols = m.ncols;
+            width = m.width;
+            rowstride = m.rowstride;
+            type = m.type;
+            data = m.data;
+            m.type = moved;
+            m.data = nullptr;
+            return *this;
+        }
+
+        ~MatGF2()
+        {
+            if (type == root)
+            {
+                CUDA_CHECK(cudaFree(data));
+            }
+        }
+
+        inline base_t *at_base(size_t r, size_t w) const
+        {
+            return data + r * rowstride + w;
+        }
+
+        void randomize(uint_fast32_t seed)
+        {
+            assert(type == root);
+            static default_random_engine e(seed);
+            static uniform_int_distribution<base_t> d;
+            base_t lastmask = base_fullmask >> (width * base_len - ncols * gf2_len);
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    *at_base(r, w) = d(e);
+                }
+                *at_base(r, width - 1) &= lastmask;
+            }
+        }
+
+        // 生成随机最简化行阶梯矩阵 前rank_col中选择nrows个主元列
+        void randomize(size_t rank_col, uint_fast32_t seed)
+        {
+            assert(nrows <= rank_col && rank_col <= ncols);
+            randomize(seed);
+            vector<size_t> pivot(rank_col);
+            iota(pivot.begin(), pivot.end(), 0);
+            random_shuffle(pivot.begin(), pivot.end());
+            pivot.resize(nrows);
+            sort(pivot.begin(), pivot.end());
+
+            vector<base_t> pivotmask(width, base_fullmask);
+            for (size_t r = 0; r < nrows; r++)
+            {
+                del(pivotmask[pivot[r] / gf2_num], pivot[r] % gf2_num);
+            }
+
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < pivot[r] / gf2_num; w++)
+                {
+                    *at_base(r, w) = base_zero;
+                }
+                base_t *now = at_base(r, pivot[r] / gf2_num);
+                *now = concat(base_zero, pivot[r] % gf2_num + 1, *now & pivotmask[pivot[r] / gf2_num]);
+                set(*now, pivot[r] % gf2_num);
+                for (size_t w = pivot[r] / gf2_num + 1; w < rank_col / gf2_num + 1; w++)
+                {
+                    *at_base(r, w) &= pivotmask[w];
+                }
+            }
+        }
+
+        bool operator==(const MatGF2 &m) const
+        {
+            if (nrows != m.nrows || ncols != m.ncols)
+            {
+                return false;
+            }
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    if (*at_base(r, w) != *m.at_base(r, w))
+                    {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+
+        bool operator==(const base_t base) const
+        {
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    if (*at_base(r, w) != base)
+                    {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+        void operator^=(const MatGF2 &m)
+        {
+            assert(nrows == m.nrows && ncols == m.ncols);
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    *at_base(r, w) ^= *m.at_base(r, w);
+                }
+            }
+        }
+        MatGF2 operator^(const MatGF2 &m) const
+        {
+            MatGF2 temp(*this);
+            temp ^= m;
+            return temp;
+        }
+
+        void gpu_addmul(const MatGF2 &a, const MatGF2 &b);
+        friend MatGF2 gpu_mul(const MatGF2 &a, const MatGF2 &b);
+
+        MatGF2 operator*(const MatGF2 &m) const
+        {
+            return gpu_mul(*this, m);
+        }
+
+        ElimResult gpu_elim();
+
+        friend ostream &operator<<(ostream &out, const MatGF2 &m);
+
+        size_t nrows, ncols, width;
+
+    private:
+        MatGF2() : nrows(0), ncols(0), width(0), rowstride(0), type(moved), data(nullptr) {}
+
+        void cpu_swap_row(size_t r1, size_t r2)
+        {
+            if (r1 == r2)
+            {
+                return;
+            }
+            base_t *p1 = at_base(r1, 0);
+            base_t *p2 = at_base(r2, 0);
+            for (size_t i = 0; i < width; i++)
+            {
+                base_t temp = p1[i];
+                p1[i] = p2[i];
+                p2[i] = temp;
+            }
+        }
+
+        size_t rowstride;
+        MatType type;
+        base_t *data;
+    };
+
+    ostream &operator<<(ostream &out, const MatGF2 &m)
+    {
+        for (size_t r = 0; r < m.nrows; r++)
+        {
+            for (size_t w = 0; w < m.width; w++)
+            {
+                bitset<gf2_num> temp(rev(*m.at_base(r, w)));
+                out << temp << " ";
+            }
+            out << endl;
+        }
+        return out;
+    }
+}
+
+#endif

include/gf2/gf2_mul.cuh

@@ -0,0 +1,83 @@
+#ifndef GF2_MUL_CUH
+#define GF2_MUL_CUH
+
+#include "gf2_mat.cuh"
+
+namespace gf2
+{
+    __global__ void gpu_addmul_kernel(base_t *a, size_t a_rowstride, base_t *tb, size_t tb_rowstride, base_t *c, size_t c_rowstride, size_t ncols, size_t width, size_t nrows)
+    {
+        size_t w = blockIdx.x * blockDim.x + threadIdx.x;
+        size_t r = blockIdx.y * blockDim.y + threadIdx.y;
+
+        if (w >= width || r >= nrows)
+        {
+            return;
+        }
+
+        base_t val = *at_base(a, a_rowstride, r, 0);
+        base_t temp = base_zero;
+        for (size_t i = 0; i < gf2_table_num; i++)
+        {
+            temp ^= *at_base(tb, tb_rowstride, i * (1 << gf2_table_len) + (val & gf2_table_mask), w);
+            val >>= gf2_table_len;
+        }
+        *at_base(c, c_rowstride, r, w) ^= temp;
+    }
+    __global__ void gpu_mktb_kernel(base_t *tb, size_t tb_rowstride, base_t *b, size_t b_rowstride, size_t tb_width)
+    {
+        size_t w = blockIdx.x * blockDim.x + threadIdx.x;
+        size_t r = blockIdx.y * blockDim.y + threadIdx.y;
+
+        if (w >= tb_width)
+        {
+            return;
+        }
+
+        base_t val = base_zero;
+        base_t idx = r & gf2_table_mask;
+        base_t st_row = (r >> gf2_table_len) * gf2_table_len;
+
+        for (size_t i = 0; i < gf2_table_len; i++)
+        {
+            if (get(idx, i) != 0)
+            {
+                val ^= *at_base(b, b_rowstride, st_row + i, w);
+            }
+        }
+        *at_base(tb, tb_rowstride, r, w) = val;
+    }
+
+    __host__ void MatGF2::gpu_addmul(const MatGF2 &a, const MatGF2 &b)
+    {
+        assert(a.ncols == b.nrows && a.nrows == nrows && b.ncols == ncols);
+        MatGF2 tb(gf2_table_num * (1 << gf2_table_len), b.ncols);
+
+        progress::ProgressBar pb("GPU MULTIPLY", a.width);
+        for (size_t w = 0; w < a.width; w++, pb.tick_display())
+        {
+            size_t size = min(base_len, a.ncols - w * base_len);
+            size_t tb_nrows = (size / gf2_table_len) * (1 << gf2_table_len) + (size % gf2_table_len == 0 ? 0 : 1 << (size % gf2_table_len));
+
+            dim3 block_tb(THREAD_X, THREAD_Y);
+            dim3 grid_tb((b.width - 1) / block_tb.x + 1, (tb_nrows - 1) / block_tb.y + 1);
+            gpu_mktb_kernel<<<grid_tb, block_tb>>>(tb.data, tb.rowstride, b.at_base(w * base_len, 0), b.rowstride, tb.width);
+            cudaDeviceSynchronize();
+
+            dim3 block(THREAD_X, THREAD_Y);
+            dim3 grid((b.width - 1) / block.x + 1, (nrows - 1) / block.y + 1);
+            gpu_addmul_kernel<<<grid, block>>>(a.at_base(0, w), a.rowstride, tb.data, tb.rowstride, data, rowstride, size, width, nrows);
+            cudaDeviceSynchronize();
+        }
+    }
+
+    __host__ MatGF2 gpu_mul(const MatGF2 &a, const MatGF2 &b)
+    {
+        assert(a.ncols == b.nrows);
+        MatGF2 c(a.nrows, b.ncols);
+        c.gpu_addmul(a, b);
+        return c;
+    }
+}
+
+#endif

include/gf256/gf256_elim.cuh

@@ -0,0 +1,179 @@
+#ifndef GF256_ELIM_CUH
+#define GF256_ELIM_CUH
+
+#include "gf256_mat.cuh"
+
+namespace gf256
+{
+    size_t cpu_elim_base(base_t *base_col, base_t base_col_len, size_t st_r, size_t w, vector<size_t> &p_col, vector<size_t> &p_row, const GF256 &gf)
+    {
+        size_t rank = 0;
+        size_t pivot[gf256_num];
+        size_t next[gf256_num];
+        for (size_t pivot_col = 0; pivot_col < gf256_num; pivot_col++)
+        {
+            for (size_t r = rank; r < base_col_len; r++)
+            {
+                for (size_t i = 0; i < rank; i++)
+                {
+                    if (next[i] == r)
+                    {
+                        base_col[r] ^= gf.mul_base(get8(base_col[r], pivot[i]), base_col[i], pivot[i] + 1);
+                        next[i]++;
+                    }
+                }
+
+                if (get8(base_col[r], pivot_col) != 0)
+                {
+                    p_col.push_back(w * gf256_num + pivot_col);
+                    p_row.push_back(st_r + r);
+                    if (r != rank)
+                    {
+                        base_t temp = base_col[rank];
+                        base_col[rank] = base_col[r];
+                        base_col[r] = temp;
+                    }
+                    base_col[rank] = concat8(base_col[rank], pivot_col + 1, gf.mul_base(gf.inv(get8(base_col[rank], pivot_col)), base_col[rank], pivot_col + 1));
+                    pivot[rank] = pivot_col;
+                    next[rank] = rank + 1;
+                    rank++;
+                    break;
+                }
+            }
+        }
+        return rank;
+    }
+
+    __global__ void gpu_mksrc_kernel(base_t *src, size_t s_rowstride, base_t *spL, size_t src_rank, size_t width)
+    {
+        size_t w = blockIdx.x * blockDim.x + threadIdx.x;
+        if (w >= width)
+        {
+            return;
+        }
+        base_t temp[gf256_num];
+        for (size_t r = 0; r < src_rank; r++)
+        {
+            temp[r] = *at_base(src, s_rowstride, r, w);
+        }
+        for (size_t r = 0; r < src_rank; r++)
+        {
+            for (size_t i = 0; i < r; i++)
+            {
+                temp[r] ^= mul_base(get8(spL[r], i), temp[i]);
+            }
+            temp[r] = mul_base(get8(spL[r], r), temp[r]);
+        }
+        for (size_t rr = 1; rr < src_rank; rr++)
+        {
+            size_t r = src_rank - 1 - rr;
+            for (size_t i = r + 1; i < src_rank; i++)
+            {
+                temp[r] ^= mul_base(get8(spL[r], i), temp[i]);
+            }
+        }
+        for (size_t r = 0; r < src_rank; r++)
+        {
+            *at_base(src, s_rowstride, r, w) = temp[r];
+        }
+    }
+
+    __global__ void gpu_elim_kernel(base_t *idx, base_t *tb, size_t tb_rowstride, base_t *data, size_t rowstride, size_t rank, base_t pivot_base, size_t st_skip, size_t width, size_t nrows)
+    {
+        size_t w = blockIdx.x * blockDim.x + threadIdx.x;
+        size_t r = blockIdx.y * blockDim.y + threadIdx.y;
+
+        if (w >= width || r >= nrows || (r >= st_skip && r < st_skip + rank))
+        {
+            return;
+        }
+
+        base_t val = idx[r];
+        base_t temp = base_zero;
+        for (size_t i = 0; i < rank; i++)
+        {
+            temp ^= *at_base(tb, tb_rowstride, i * (1 << gf256_len) + get8(val, get8(pivot_base, i)), w);
+        }
+        *at_base(data, rowstride, r, w) ^= temp;
+    }
+
+    __managed__ base_t spL[gf256_num];
+
+    __host__ ElimResult MatGF256::gpu_elim(const GF256 &gf)
+    {
+        gf.cpy_to_constant();
+        MatGF256 tb(gf256_num * (1 << gf256_len), ncols);
+
+        base_t *base_col;
+        cudaMallocManaged(&base_col, nrows * sizeof(base_t));
+        base_t *idx;
+        cudaMallocManaged(&idx, nrows * sizeof(base_t));
+
+        size_t rank = 0;
+        vector<size_t> p_col, p_row;
+
+        progress::ProgressBar pb("GPU ELIMINATE", width);
+        for (size_t w = 0; w < width; w++, pb.tick_display())
+        {
+            CUDA_CHECK(cudaMemcpy2D(base_col + rank, sizeof(base_t), at_base(rank, w), rowstride * sizeof(base_t), sizeof(base_t), nrows - rank, cudaMemcpyDefault));
+
+            size_t src_rank = cpu_elim_base(base_col + rank, nrows - rank, rank, w, p_col, p_row, gf);
+
+            if (src_rank == 0)
+            {
+                continue;
+            }
+
+            for (size_t i = 0; i < src_rank; i++)
+            {
+                cpu_swap_row(rank + i, p_row[rank + i]);
+                spL[i] = base_zero;
+            }
+
+            base_t pivot_base = base_zero;
+            for (size_t r = 0; r < src_rank; r++)
+            {
+                size_t loc = (p_col[rank + r] - w * gf256_num);
+                set8(spL[r], gf.inv(get8(base_col[rank + r], loc)), r);
+                for (size_t i = 0; i < r; i++)
+                {
+                    set8(spL[i], get8(base_col[rank + i], loc), r);
+                }
+                for (size_t i = r + 1; i < src_rank; i++)
+                {
+                    set8(spL[i], get8(base_col[rank + i], loc), r);
+                }
+                set8(pivot_base, loc, r);
+            }
+
+            dim3 block_src(THREAD_X);
+            dim3 grid_src((width - w - 1) / block_src.x + 1);
+            gpu_mksrc_kernel<<<grid_src, block_src>>>(at_base(rank, w), rowstride, spL, src_rank, width);
+            cudaDeviceSynchronize();
+
+            dim3 block_tb(THREAD_X, THREAD_Y);
+            dim3 grid_tb((width - w - 1) / block_tb.x + 1, (src_rank * (1 << gf256_len) - 1) / block_tb.y + 1);
+            gpu_mktb_kernel<<<grid_tb, block_tb>>>(tb.data, tb.rowstride, at_base(rank, w), rowstride, tb.width);
+            cudaDeviceSynchronize();
+
+            CUDA_CHECK(cudaMemcpy2D(idx, sizeof(base_t), at_base(0, w), rowstride * sizeof(base_t), sizeof(base_t), nrows, cudaMemcpyDefault));
+
+            dim3 block(THREAD_X, THREAD_Y);
+            dim3 grid((width - w - 1) / block.x + 1, (nrows - 1) / block.y + 1);
+            gpu_elim_kernel<<<grid, block>>>(idx, tb.data, tb.rowstride, at_base(0, w), rowstride, src_rank, pivot_base, rank, width - w, nrows);
+            cudaDeviceSynchronize();
+
+            rank += src_rank;
+
+            if (rank == nrows)
+            {
+                break;
+            }
+        }
+        cudaFree(base_col);
+        cudaFree(idx);
+        return {rank, p_col, p_row};
+    }
+}
+
+#endif

include/gf256/gf256_header.cuh

@@ -0,0 +1,199 @@
+#ifndef GF256_HEADER_CUH
+#define GF256_HEADER_CUH
+
+#include "../header.cuh"
+#include <set>
+
+namespace gf256
+{
+    using gf256_t = uint8_t;
+
+    static const size_t gf256_len = sizeof(gf256_t) * 8;
+    static const size_t gf256_num = base_len / gf256_len;
+
+    static const gf256_t gf256_zero = (gf256_t)0x00;
+    static const gf256_t gf256_one = (gf256_t)0x01;
+    static const gf256_t gf256_fullmask = (gf256_t)0xFF;
+
+    static const base_t gf256_mask[8] = {
+        (base_t)0x00'00'00'00'00'00'00'FF,
+        (base_t)0x00'00'00'00'00'00'FF'00,
+        (base_t)0x00'00'00'00'00'FF'00'00,
+        (base_t)0x00'00'00'00'FF'00'00'00,
+        (base_t)0x00'00'00'FF'00'00'00'00,
+        (base_t)0x00'00'FF'00'00'00'00'00,
+        (base_t)0x00'FF'00'00'00'00'00'00,
+        (base_t)0xFF'00'00'00'00'00'00'00};
+
+    __host__ __device__ inline size_t offset8(size_t idx)
+    {
+        return idx << 3;
+    }
+
+    __host__ __device__ inline gf256_t get8(base_t src, size_t idx)
+    {
+        return (gf256_t)(src >> offset8(idx));
+    }
+
+    // 确保set8对应位置的值为0
+    __host__ __device__ inline void set8(base_t &dst, gf256_t src, size_t idx)
+    {
+        dst |= (base_t)src << offset8(idx);
+    }
+
+    __host__ inline void del8(base_t &dst, size_t idx)
+    {
+        dst &= ~gf256_mask[idx];
+    }
+
+    __host__ inline base_t concat8(base_t dst_l, size_t idx_l, base_t dst_r)
+    {
+        if (idx_l == 0)
+        {
+            return dst_r;
+        }
+        if (idx_l == gf256_num)
+        {
+            return dst_l;
+        }
+        return (dst_l & (base_fullmask >> (base_len - offset8(idx_l)))) | (dst_r & (base_fullmask << offset8(idx_l)));
+    }
+
+    __host__ inline base_t rev8(base_t n)
+    {
+        n = (n & (base_t)0xFF'00'FF'00'FF'00'FF'00) >> 8 | (n & (base_t)0x00'FF'00'FF'00'FF'00'FF) << 8;
+        n = (n & (base_t)0xFF'FF'00'00'FF'FF'00'00) >> 16 | (n & (base_t)0x00'00'FF'FF'00'00'FF'FF) << 16;
+        return n >> 32 | n << 32;
+    }
+
+    __constant__ gf256_t d_mul_table[1 << gf256_len][1 << gf256_len];
+
+    __device__ inline base_t mul_base(const gf256_t val, const base_t base)
+    {
+        if (val == 0)
+        {
+            return base_zero;
+        }
+        base_t temp = base_zero;
+        for (size_t i = 0; i < gf256_len; i++)
+        {
+            set8(temp, d_mul_table[val][get8(base, i)], i);
+        }
+        return temp;
+    }
+
+    __global__ void gpu_mktb_kernel(base_t *tb, size_t tb_rowstride, base_t *src, size_t s_rowstride, size_t width)
+    {
+        size_t w = blockIdx.x * blockDim.x + threadIdx.x;
+        size_t r = blockIdx.y * blockDim.y + threadIdx.y;
+
+        if (w >= width)
+        {
+            return;
+        }
+
+        gf256_t val = get8(r, 0);
+        base_t s = *at_base(src, s_rowstride, get8(r, 1), w);
+        base_t d = mul_base(val, s);
+        *at_base(tb, tb_rowstride, r, w) = d;
+    }
+
+    static const set<base_t> irreducible_polynomials_degree_08{0x11b, 0x11d, 0x12b, 0x12d, 0x139, 0x13f, 0x14d, 0x15f, 0x163, 0x165, 0x169, 0x171, 0x177, 0x17b, 0x187, 0x18b, 0x18d, 0x19f, 0x1a3, 0x1a9, 0x1b1, 0x1bd, 0x1c3, 0x1cf, 0x1d7, 0x1dd, 0x1e7, 0x1f3, 0x1f5, 0x1f9};
+
+    class GF256
+    {
+    public:
+        GF256(base_t poly)
+        {
+            assert(irreducible_polynomials_degree_08.count(poly) == 1);
+            this->polynomial = poly;
+            for (size_t x = 0; x < (1 << gf256_len); x++)
+            {
+                mul_table[x][gf256_zero] = gf256_zero;
+                for (size_t d = 0; d < gf256_len; d++)
+                {
+                    gf256_t val = shift_left(x, d);
+                    for (size_t y = (1 << d); y < (1 << (d + 1)); y++)
+                    {
+                        mul_table[x][y] = val ^ mul_table[x][y ^ (1 << d)];
+                        if (mul_table[x][y] == gf256_one)
+                        {
+                            inv_table[x] = y;
+                        }
+                    }
+                }
+            }
+            inv_table[gf256_zero] = gf256_zero;
+        }
+
+        gf256_t mul(const gf256_t x, const gf256_t y) const
+        {
+            return mul_table[x][y];
+        }
+
+        base_t mul_base(const gf256_t val, const base_t base, const size_t offset = 0) const
+        {
+            base_t temp = base_zero;
+            for (size_t i = offset; i < gf256_num; i++)
+            {
+                set8(temp, mul(val, get8(base, i)), i);
+            }
+            return temp;
+        }
+
+        gf256_t inv(gf256_t x) const
+        {
+            return inv_table[x];
+        }
+
+        base_t poly(void) const
+        {
+            return polynomial;
+        }
+
+        inline cudaError_t cpy_to_constant() const
+        {
+            return cudaMemcpyToSymbol(d_mul_table, mul_table, (1 << gf256_len) * (1 << gf256_len) * sizeof(gf256_t));
+        }
+
+        friend ostream &operator<<(ostream &out, const GF256 &gf);
+
+        GF256() = delete;
+        GF256(const GF256 &) = delete;
+        GF256(GF256 &&) = delete;
+        GF256 &operator=(const GF256 &) = delete;
+        GF256 &operator=(GF256 &&) = delete;
+
+    private:
+        gf256_t shift_left(gf256_t x, size_t d)
+        {
+            base_t temp = (base_t)x << d;
+            for (size_t i = gf256_len - 1 + d; i > gf256_len - 1; i--)
+            {
+                if (temp & (1 << i))
+                {
+                    temp ^= polynomial << (i - gf256_len);
+                }
+            }
+            return temp;
+        }
+
+        base_t polynomial;
+        gf256_t inv_table[1 << gf256_num];
+        gf256_t mul_table[1 << gf256_num][1 << gf256_num];
+    };
+
+    ostream &operator<<(ostream &out, const GF256 &gf)
+    {
+        for (size_t x = 0; x < 1 << gf256_len; x++)
+        {
+            for (size_t y = 0; y < 1 << gf256_len; y++)
+            {
+                printf("%02X ", gf.mul_table[x][y]);
+            }
+            printf("\n");
+        }
+        return out;
+    }
+}
+#endif

include/gf256/gf256_mat.cuh

@@ -0,0 +1,242 @@
+#ifndef GF256_MAT_CUH
+#define GF256_MAT_CUH
+
+#include "gf256_header.cuh"
+
+#include <random>
+#include <vector>
+#include <algorithm>
+
+namespace gf256
+{
+    struct ElimResult
+    {
+        size_t rank;
+        vector<size_t> pivot;
+        vector<size_t> swap_row;
+    };
+
+    class MatGF256
+    {
+    public:
+        enum MatType
+        {
+            root,
+            window,
+            moved,
+        };
+        // 只能构造root矩阵
+        MatGF256(size_t nrows, size_t ncols) : nrows(nrows), ncols(ncols), type(root)
+        {
+            width = (ncols - 1) / gf256_num + 1;
+            rowstride = ((width - 1) / 4 + 1) * 4; // 以32字节（4*64bit）为单位对齐
+            CUDA_CHECK(cudaMallocManaged((void **)&data, nrows * rowstride * sizeof(base_t)));
+            CUDA_CHECK(cudaMemset(data, 0, nrows * rowstride * sizeof(base_t)));
+        }
+        // 只能以base_t为单位建立window矩阵
+        MatGF256(const MatGF256 &src, size_t begin_ri, size_t begin_wi, size_t end_rj, size_t end_wj) : nrows(end_rj - begin_ri), ncols((end_wj == src.width ? src.ncols : end_wj * gf256_num) - begin_wi * gf256_num), width(end_wj - begin_wi), rowstride(src.rowstride), type(window), data(src.at_base(begin_ri, begin_wi))
+        {
+            assert(begin_ri < end_rj && end_rj <= src.nrows && begin_wi < end_wj && end_wj <= src.width);
+        }
+        // 只能拷贝构造root矩阵
+        MatGF256(const MatGF256 &m) : MatGF256(m.nrows, m.ncols)
+        {
+            CUDA_CHECK(cudaMemcpy2D(data, rowstride * sizeof(base_t), m.data, m.rowstride * sizeof(base_t), m.width * sizeof(base_t), nrows, cudaMemcpyDefault));
+        }
+        MatGF256(MatGF256 &&m) noexcept : nrows(m.nrows), ncols(m.ncols), width(m.width), rowstride(m.rowstride), type(m.type), data(m.data)
+        {
+            m.type = moved;
+            m.data = nullptr;
+        }
+        MatGF256 &operator=(const MatGF256 &m)
+        {
+            if (this == &m)
+            {
+                return *this;
+            }
+            assert(nrows == m.nrows && ncols == m.ncols);
+            CUDA_CHECK(cudaMemcpy2D(data, rowstride * sizeof(base_t), m.data, m.rowstride * sizeof(base_t), m.width * sizeof(base_t), nrows, cudaMemcpyDefault));
+            return *this;
+        }
+        MatGF256 &operator=(MatGF256 &&m) noexcept
+        {
+            if (this == &m)
+            {
+                return *this;
+            }
+            if (type == root)
+            {
+                CUDA_CHECK(cudaFree(data));
+            }
+            nrows = m.nrows;
+            ncols = m.ncols;
+            width = m.width;
+            rowstride = m.rowstride;
+            type = m.type;
+            data = m.data;
+            m.type = moved;
+            m.data = nullptr;
+            return *this;
+        }
+
+        ~MatGF256()
+        {
+            if (type == root)
+            {
+                CUDA_CHECK(cudaFree(data));
+            }
+        }
+
+        inline base_t *at_base(size_t r, size_t w) const
+        {
+            return data + r * rowstride + w;
+        }
+
+        void randomize(uint_fast32_t seed)
+        {
+            assert(type == root);
+            static default_random_engine e(seed);
+            static uniform_int_distribution<base_t> d;
+            base_t lastmask = base_fullmask >> (width * base_len - ncols * gf256_len);
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    *at_base(r, w) = d(e);
+                }
+                *at_base(r, width - 1) &= lastmask;
+            }
+        }
+
+        // 生成随机最简化行阶梯矩阵 前rank_col中选择nrows个主元列
+        void randomize(size_t rank_col, uint_fast32_t seed)
+        {
+            assert(nrows <= rank_col && rank_col <= ncols);
+            randomize(seed);
+            vector<size_t> pivot(rank_col);
+            iota(pivot.begin(), pivot.end(), 0);
+            random_shuffle(pivot.begin(), pivot.end());
+            pivot.resize(nrows);
+            sort(pivot.begin(), pivot.end());
+
+            vector<base_t> pivotmask(width, base_fullmask);
+            for (size_t r = 0; r < nrows; r++)
+            {
+                del8(pivotmask[pivot[r] / gf256_num], pivot[r] % gf256_num);
+            }
+
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < pivot[r] / gf256_num; w++)
+                {
+                    *at_base(r, w) = base_zero;
+                }
+                base_t *now = at_base(r, pivot[r] / gf256_num);
+                *now = concat8(base_zero, pivot[r] % gf256_num + 1, *now & pivotmask[pivot[r] / gf256_num]);
+                set8(*now, gf256_one, pivot[r] % gf256_num);
+                for (size_t w = pivot[r] / gf256_num + 1; w < rank_col / gf256_num + 1; w++)
+                {
+                    *at_base(r, w) &= pivotmask[w];
+                }
+            }
+        }
+
+        bool operator==(const MatGF256 &m) const
+        {
+            if (nrows != m.nrows || ncols != m.ncols)
+            {
+                return false;
+            }
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    if (*at_base(r, w) != *m.at_base(r, w))
+                    {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+
+        bool operator==(const base_t base) const
+        {
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    if (*at_base(r, w) != base)
+                    {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+        void operator^=(const MatGF256 &m)
+        {
+            assert(nrows == m.nrows && ncols == m.ncols);
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    *at_base(r, w) ^= *m.at_base(r, w);
+                }
+            }
+        }
+        MatGF256 operator^(const MatGF256 &m) const
+        {
+            MatGF256 temp(*this);
+            temp ^= m;
+            return temp;
+        }
+
+        void gpu_addmul(const MatGF256 &a, const MatGF256 &b, const GF256 &gf);
+        friend MatGF256 gpu_mul(const MatGF256 &a, const MatGF256 &b, const GF256 &gf);
+
+        ElimResult gpu_elim(const GF256 &gf);
+
+        friend ostream &operator<<(ostream &out, const MatGF256 &m);
+
+        size_t nrows, ncols, width;
+
+    private:
+        MatGF256() : nrows(0), ncols(0), width(0), rowstride(0), type(moved), data(nullptr) {}
+
+        void cpu_swap_row(size_t r1, size_t r2)
+        {
+            if (r1 == r2)
+            {
+                return;
+            }
+            base_t *p1 = at_base(r1, 0);
+            base_t *p2 = at_base(r2, 0);
+            for (size_t i = 0; i < width; i++)
+            {
+                base_t temp = p1[i];
+                p1[i] = p2[i];
+                p2[i] = temp;
+            }
+        }
+
+        size_t rowstride;
+        MatType type;
+        base_t *data;
+    };
+
+    ostream &operator<<(ostream &out, const MatGF256 &m)
+    {
+        for (size_t r = 0; r < m.nrows; r++)
+        {
+            for (size_t w = 0; w < m.width; w++)
+            {
+                printf("%016lX ", rev8(*m.at_base(r, w)));
+            }
+            printf("\n");
+        }
+        return out;
+    }
+}
+
+#endif

include/gf256/gf256_mul.cuh

@@ -0,0 +1,59 @@
+#ifndef GF256_MUL_CUH
+#define GF256_MUL_CUH
+
+#include "gf256_mat.cuh"
+
+namespace gf256
+{
+    __global__ void gpu_addmul_kernel(base_t *a, size_t a_rowstride, base_t *tb, size_t tb_rowstride, base_t *c, size_t c_rowstride, size_t tb_num, size_t width, size_t nrows)
+    {
+        size_t w = blockIdx.x * blockDim.x + threadIdx.x;
+        size_t r = blockIdx.y * blockDim.y + threadIdx.y;
+
+        if (w >= width || r >= nrows)
+        {
+            return;
+        }
+
+        base_t val = *at_base(a, a_rowstride, r, 0);
+        base_t temp = base_zero;
+        for (size_t i = 0; i < tb_num; i++)
+        {
+            temp ^= *at_base(tb, tb_rowstride, i * (1 << gf256_len) + get8(val, i), w);
+        }
+        *at_base(c, c_rowstride, r, w) ^= temp;
+    }
+
+    __host__ void MatGF256::gpu_addmul(const MatGF256 &a, const MatGF256 &b, const GF256 &gf)
+    {
+        assert(a.ncols == b.nrows && a.nrows == nrows && b.ncols == ncols);
+        gf.cpy_to_constant();
+        MatGF256 tb(gf256_num * (1 << gf256_len), b.ncols);
+
+        progress::ProgressBar pb("GPU MULTIPLY", a.width);
+        for (size_t w = 0; w < a.width; w++, pb.tick_display())
+        {
+            size_t tb_num = min(gf256_num, a.ncols - w * gf256_num);
+
+            dim3 block_tb(THREAD_X, THREAD_Y);
+            dim3 grid_tb((b.width - 1) / block_tb.x + 1, (tb_num * (1 << gf256_len) - 1) / block_tb.y + 1);
+            gpu_mktb_kernel<<<grid_tb, block_tb>>>(tb.data, tb.rowstride, b.at_base(w * gf256_num, 0), b.rowstride, tb.width);
+            cudaDeviceSynchronize();
+
+            dim3 block(THREAD_X, THREAD_Y);
+            dim3 grid((b.width - 1) / block.x + 1, (nrows - 1) / block.y + 1);
+            gpu_addmul_kernel<<<grid, block>>>(a.at_base(0, w), a.rowstride, tb.data, tb.rowstride, data, rowstride, tb_num, width, nrows);
+            cudaDeviceSynchronize();
+        }
+    }
+
+    __host__ MatGF256 gpu_mul(const MatGF256 &a, const MatGF256 &b, const GF256 &gf)
+    {
+        assert(a.ncols == b.nrows);
+        MatGF256 c(a.nrows, b.ncols);
+        c.gpu_addmul(a, b, gf);
+        return c;
+    }
+}
+
+#endif

include/gf28.cuh

@@ -1,95 +0,0 @@
-#ifndef GF28_CUH
-#define GF28_CUH
-
-#include "header.cuh"
-
-static const set<base_t> irreducible_polynomials_degree_08{0x11b, 0x11d, 0x12b, 0x12d, 0x139, 0x13f, 0x14d, 0x15f, 0x163, 0x165, 0x169, 0x171, 0x177, 0x17b, 0x187, 0x18b, 0x18d, 0x19f, 0x1a3, 0x1a9, 0x1b1, 0x1bd, 0x1c3, 0x1cf, 0x1d7, 0x1dd, 0x1e7, 0x1f3, 0x1f5, 0x1f9};
-
-class GF28
-{
-public:
-    GF28(base_t poly)
-    {
-        assert(irreducible_polynomials_degree_08.count(poly) == 1);
-        this->poly = poly;
-        for (size_t x = 0; x < (1 << base_deg); x++)
-        {
-            mul_table[x][gf28_zero] = gf28_zero;
-            for (size_t d = 0; d < base_deg; d++)
-            {
-                gf28_t val = shift_left(x, d);
-                for (size_t y = (1 << d); y < (1 << (d + 1)); y++)
-                {
-                    mul_table[x][y] = val ^ mul_table[x][y ^ (1 << d)];
-                    if (mul_table[x][y] == gf28_one)
-                    {
-                        inv_table[x] = y;
-                    }
-                }
-            }
-        }
-        inv_table[gf28_zero] = gf28_zero;
-    }
-
-    gf28_t mul(const gf28_t x, const gf28_t y) const
-    {
-        return mul_table[x][y];
-    }
-
-    base_t mul_base(const gf28_t val, const base_t base, const size_t offset = 0) const
-    {
-        base_t temp = base_zero;
-        for (size_t i = offset; i < base_num; i++)
-        {
-            set8(temp, mul(val, get8(base, i)), i);
-        }
-        return temp;
-    }
-
-    gf28_t inv(gf28_t x)
-    {
-        return inv_table[x];
-    }
-
-    friend ostream &operator<<(ostream &out, const GF28 &gf);
-
-    GF28() = delete;
-    GF28(const GF28 &) = delete;
-    GF28(GF28 &&) = delete;
-    GF28 &operator=(const GF28 &) = delete;
-    GF28 &operator=(GF28 &&) = delete;
-
-    gf28_t mul_table[1 << base_deg][1 << base_deg];
-
-private:
-    gf28_t shift_left(gf28_t x, size_t d)
-    {
-        base_t temp = (base_t)x << d;
-        for (size_t i = base_deg - 1 + d; i > base_deg - 1; i--)
-        {
-            if (temp & (1 << i))
-            {
-                temp ^= poly << (i - base_deg);
-            }
-        }
-        return temp;
-    }
-
-    base_t poly;
-    gf28_t inv_table[1 << base_deg];
-};
-
-ostream &operator<<(ostream &out, const GF28 &gf)
-{
-    for (size_t x = 0; x < 1 << base_deg; x++)
-    {
-        for (size_t y = 0; y < 1 << base_deg; y++)
-        {
-            printf("%02X ", gf.mul_table[x][y]);
-        }
-        printf("\n");
-    }
-    return out;
-}
-
-#endif

include/gfp/gfp_elim.cuh

@@ -0,0 +1,218 @@
+#ifndef GFP_ELIM_CUH
+#define GFP_ELIM_CUH
+
+#include "gfp_mat.cuh"
+
+static const size_t gfp_block_size = 8;
+
+namespace gfp
+{
+    ElimResult MatGFP::cpu_elim()
+    {
+        init_inv_table();
+        size_t rank = 0;
+        vector<size_t> pivot;
+        vector<size_t> swap_row;
+        for (size_t w = 0; w < width; w++)
+        {
+            if (*at_base(rank, w) != gfp_zero)
+            {
+                swap_row.push_back(rank);
+            }
+            else
+            {
+                for (size_t r = rank + 1; r < nrows; r++)
+                {
+                    if (*at_base(r, w) != gfp_zero)
+                    {
+                        swap_row.push_back(r);
+                        cpu_swap_row(rank, r);
+                        break;
+                    }
+                }
+                continue;
+            }
+            pivot.push_back(w);
+
+            gfp_t inv = gfp_inv(*at_base(rank, w));
+
+            for (size_t wi = w; wi < width; wi++)
+            {
+                *at_base(rank, wi) = (*at_base(rank, wi) * inv) % gfprime;
+            }
+
+            for (size_t r = 0; r < nrows; r++)
+            {
+                if (r == rank)
+                {
+                    continue;
+                }
+
+                gfp_t mul = *at_base(r, w);
+
+                for (size_t wi = w; wi < width; wi++)
+                {
+                    *at_base(r, wi) = (*at_base(r, wi) + (*at_base(rank, wi) * (gfprime - mul))) % gfprime;
+                }
+            }
+
+            rank++;
+            if (rank == nrows)
+            {
+                break;
+            }
+        }
+        return {rank, pivot, swap_row};
+    }
+
+    size_t cpu_elim_base(MatGFP &submat, base_t submat_nrows, size_t st_r, size_t w, vector<size_t> &p_col, vector<size_t> &p_row)
+    {
+        size_t rank = 0;
+        size_t pivot[gfp_block_size];
+        size_t next[gfp_block_size];
+        for (size_t pivot_col = 0; pivot_col < gfp_block_size; pivot_col++)
+        {
+            for (size_t r = rank; r < submat_nrows; r++)
+            {
+                for (size_t i = 0; i < rank; i++)
+                {
+                    if (next[i] == r)
+                    {
+                        gfp_t mul = *submat.at_base(r, pivot[i]);
+                        for (size_t j = pivot[i] + 1; j < gfp_block_size; j++)
+                        {
+                            *submat.at_base(r, j) = (*submat.at_base(r, j) + *submat.at_base(i, j) * (gfprime - mul)) % gfprime;
+                        }
+                        next[i]++;
+                    }
+                }
+
+                if (*submat.at_base(r, pivot_col) != 0)
+                {
+                    p_col.push_back(w + pivot_col);
+                    p_row.push_back(st_r + r);
+                    if (r != rank)
+                    {
+                        submat.cpu_swap_row(rank, r);
+                    }
+                    gfp_t inv = gfp_inv(*submat.at_base(rank, pivot_col));
+                    for (size_t j = pivot_col + 1; j < gfp_block_size; j++)
+                    {
+                        *submat.at_base(rank, j) = (inv * (*submat.at_base(rank, j))) % gfprime;
+                    }
+                    pivot[rank] = pivot_col;
+                    next[rank] = rank + 1;
+                    rank++;
+                    break;
+                }
+            }
+        }
+        return rank;
+    }
+
+    __global__ void gpu_elim_kernel(gfp_t *__restrict__ idx, const size_t idx_rs, gfp_t *__restrict__ data, const size_t rs, const size_t st_src, const size_t rank, const gfp_t pivot_base, const size_t nrows, const size_t ncols, const size_t nsteps)
+    {
+
+        const unsigned int bx = blockIdx.x;
+        const unsigned int by = blockIdx.y;
+        const unsigned int tx = threadIdx.x;
+        const unsigned int ty = threadIdx.y;
+        const unsigned int tid = ty * blockDim.x + tx;
+
+        __shared__ gfp_t s_idx[StepSize][BlockRow];
+        __shared__ gfp_t s_src[StepSize][BlockCol];
+
+        gfp_t tmp_c[BlockRow / THREAD_Y][BlockCol / THREAD_X] = {0};
+        for (int k = tid; k < StepSize * BlockRow; k += blockDim.x * blockDim.y)
+        {
+            const int a_r = k / StepSize;
+            const int a_c = k % StepSize;
+            s_idx[a_c][a_r] = (by * BlockRow + a_r < nrows && a_c < nsteps) ? *at_base(idx, idx_rs, by * BlockRow + a_r, a_c) : 0;
+            const int b_r = k / BlockCol;
+            const int b_c = k % BlockCol;
+            s_src[b_r][b_c] = (b_r < rank && bx * BlockCol + b_c < ncols) ? *at_base(data, rs, st_src + b_r, bx * BlockCol + b_c) : 0;
+            s_src[b_r][b_c] = s_src[b_r][b_c] ? gfprime - s_src[b_r][b_c] : 0;
+        }
+        __syncthreads();
+        for (int k = 0; k < rank; k++)
+        {
+            for (int j = 0; j < BlockRow / THREAD_Y; j++)
+            {
+                for (int i = 0; i < BlockCol / THREAD_X; i++)
+                {
+                    tmp_c[j][i] += (s_idx[get4(pivot_base, k)][j * THREAD_Y + ty] * s_src[k][i * THREAD_X + tx]) % gfprime;
+                }
+            }
+        }
+
+        for (int j = 0; j < BlockRow / THREAD_Y; j++)
+        {
+            for (int i = 0; i < BlockCol / THREAD_X; i++)
+            {
+                if (by * BlockRow + j * THREAD_Y + ty < nrows && (by * BlockRow + j * THREAD_Y + ty < st_src || by * BlockRow + j * THREAD_Y + ty >= st_src + rank) && bx * BlockCol + i * THREAD_X + tx < ncols)
+                {
+                    gfp_t *dst = at_base(data, rs, by * BlockRow + j * THREAD_Y + ty, bx * BlockCol + i * THREAD_X + tx);
+                    *dst = (*dst + tmp_c[j][i]) % gfprime;
+                }
+            }
+        }
+    }
+
+    __host__ ElimResult MatGFP::gpu_elim()
+    {
+        init_inv_table();
+
+        MatGFP submat(nrows, gfp_block_size);
+
+        size_t rank = 0;
+        vector<size_t> p_col, p_row;
+
+        progress::ProgressBar pb("GPU ELIMINATE", (width - 1) / gfp_block_size + 1);
+        for (size_t w = 0; w < (width - 1) / gfp_block_size + 1; w++, pb.tick_display())
+        {
+            size_t st_w = w * gfp_block_size;
+            size_t ed_w = min(st_w + gfp_block_size, width);
+
+            submat.deepcopy(*this, rank, st_w, nrows, ed_w);
+
+            size_t src_rank = cpu_elim_base(submat, nrows - rank, rank, st_w, p_col, p_row);
+
+            if (src_rank == 0)
+            {
+                continue;
+            }
+            gfp_t pivot_base = gfp_zero;
+            for (size_t r = rank; r < rank + src_rank; r++)
+            {
+                cpu_swap_row(r, p_row[r]);
+                for (size_t i = rank; i < r; i++)
+                {
+                    cpu_addmul_row(r, i, gfprime - *at_base(r, p_col[i]), p_col[i]);
+                }
+                cpu_mul_row(r, gfp_inv(*at_base(r, p_col[r])), p_col[r]);
+                for (size_t i = rank; i < r; i++)
+                {
+                    cpu_addmul_row(i, r, gfprime - *at_base(i, p_col[r]), p_col[r]);
+                }
+                set4(pivot_base, p_col[r] - st_w, r - rank);
+            }
+
+            submat.deepcopy(*this, 0, st_w, nrows, ed_w);
+
+            dim3 block(THREAD_X, THREAD_Y);
+            dim3 grid((width - st_w - 1) / block.x + 1, (nrows - 1) / block.y + 1);
+            gpu_elim_kernel<<<grid, block>>>(submat.data, submat.rowstride, at_base(0, st_w), rowstride, rank, src_rank, pivot_base, nrows, width - st_w, gfp_block_size);
+            cudaDeviceSynchronize();
+
+            rank += src_rank;
+
+            if (rank == nrows)
+            {
+                break;
+            }
+        }
+        return {rank, p_col, p_row};
+    }
+}
+
+#endif

include/gfp/gfp_header.cuh

@@ -0,0 +1,50 @@
+#ifndef GFP_HEADER_CUH
+#define GFP_HEADER_CUH
+
+#include "../header.cuh"
+
+namespace gfp
+{
+    using gfp_t = uint32_t;
+
+    static const gfp_t gfprime = 65521;
+
+    static const gfp_t gfp_zero = (gfp_t)0;
+    static const gfp_t gfp_one = (gfp_t)1;
+    static const gfp_t gfp_fullmask = (gfp_t)0xFF'FF;
+
+    __managed__ gfp_t gfp_inv_table[gfprime];
+
+    void init_inv_table()
+    {
+        gfp_inv_table[0] = 0;
+        gfp_inv_table[1] = 1;
+        for (int i = 2; i < gfprime; ++i)
+        {
+            gfp_inv_table[i] = (gfprime - gfprime / i) * gfp_inv_table[gfprime % i] % gfprime;
+        }
+    }
+    gfp_t gfp_inv(gfp_t a)
+    {
+        return gfp_inv_table[a];
+    }
+
+    __host__ __device__ inline size_t offset4(size_t idx)
+    {
+        return idx << 2;
+    }
+
+    __host__ __device__ inline size_t get4(gfp_t src, size_t idx)
+    {
+        return (src >> offset4(idx)) & 0xF;
+    }
+
+    // 确保set8对应位置的值为0
+    __host__ __device__ inline void set4(gfp_t &dst, size_t src, size_t idx)
+    {
+        dst |= (gfp_t)src << offset4(idx);
+    }
+
+}
+
+#endif

include/gfp/gfp_mat.cuh

@@ -0,0 +1,288 @@
+#ifndef GFP_MAT_CUH
+#define GFP_MAT_CUH
+
+#include "gfp_header.cuh"
+
+#include <random>
+#include <vector>
+#include <algorithm>
+
+namespace gfp
+{
+    struct ElimResult
+    {
+        size_t rank;
+        vector<size_t> pivot;
+        vector<size_t> swap_row;
+    };
+
+    class MatGFP
+    {
+    public:
+        enum MatType
+        {
+            root,
+            window,
+            moved,
+        };
+        // 只能构造root矩阵
+        MatGFP(size_t nrows, size_t ncols) : nrows(nrows), ncols(ncols), type(root)
+        {
+            width = ncols;
+            rowstride = ((width - 1) / 8 + 1) * 8; // 以32字节（8*32bit）为单位对齐
+            CUDA_CHECK(cudaMallocManaged((void **)&data, nrows * rowstride * sizeof(gfp_t)));
+            CUDA_CHECK(cudaMemset(data, 0, nrows * rowstride * sizeof(gfp_t)));
+        }
+        // 只能以gfp_t为单位建立window矩阵
+        MatGFP(const MatGFP &src, size_t begin_ri, size_t begin_wi, size_t end_rj, size_t end_wj) : nrows(end_rj - begin_ri), ncols(end_wj - begin_wi), width(end_wj - begin_wi), rowstride(src.rowstride), type(window), data(src.at_base(begin_ri, begin_wi))
+        {
+            assert(begin_ri < end_rj && end_rj <= src.nrows && begin_wi < end_wj && end_wj <= src.width);
+        }
+        // 只能拷贝构造root矩阵
+        MatGFP(const MatGFP &m) : MatGFP(m.nrows, m.ncols)
+        {
+            CUDA_CHECK(cudaMemcpy2D(data, rowstride * sizeof(gfp_t), m.data, m.rowstride * sizeof(gfp_t), m.width * sizeof(gfp_t), nrows, cudaMemcpyDefault));
+        }
+        MatGFP(MatGFP &&m) noexcept : nrows(m.nrows), ncols(m.ncols), width(m.width), rowstride(m.rowstride), type(m.type), data(m.data)
+        {
+            m.type = moved;
+            m.data = nullptr;
+        }
+        MatGFP &operator=(const MatGFP &m)
+        {
+            if (this == &m)
+            {
+                return *this;
+            }
+            assert(nrows == m.nrows && ncols == m.ncols);
+            CUDA_CHECK(cudaMemcpy2D(data, rowstride * sizeof(gfp_t), m.data, m.rowstride * sizeof(gfp_t), m.width * sizeof(gfp_t), nrows, cudaMemcpyDefault));
+            return *this;
+        }
+        MatGFP &operator=(MatGFP &&m) noexcept
+        {
+            if (this == &m)
+            {
+                return *this;
+            }
+            if (type == root)
+            {
+                CUDA_CHECK(cudaFree(data));
+            }
+            nrows = m.nrows;
+            ncols = m.ncols;
+            width = m.width;
+            rowstride = m.rowstride;
+            type = m.type;
+            data = m.data;
+            m.type = moved;
+            m.data = nullptr;
+            return *this;
+        }
+
+        ~MatGFP()
+        {
+            if (type == root)
+            {
+                // cout << nrows << " " << ncols << endl;
+                // cout << *data << endl;
+                CUDA_CHECK(cudaFree(data));
+            }
+        }
+
+        inline gfp_t *at_base(size_t r, size_t w) const
+        {
+            return data + r * rowstride + w;
+        }
+
+        void randomize(uint_fast32_t seed)
+        {
+            assert(type == root);
+            static default_random_engine e(seed);
+            static uniform_int_distribution<gfp_t> d;
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    *at_base(r, w) = d(e) % gfprime;
+                }
+            }
+        }
+
+        void deepcopy(const MatGFP &src, size_t begin_ri, size_t begin_wi, size_t end_rj, size_t end_wj)
+        {
+            assert(end_rj - begin_ri <= nrows && end_wj - begin_wi <= width);
+            CUDA_CHECK(cudaMemcpy2D(data, rowstride * sizeof(gfp_t), src.at_base(begin_ri, begin_wi), src.rowstride * sizeof(gfp_t), (end_wj - begin_wi) * sizeof(gfp_t), end_rj - begin_ri, cudaMemcpyDefault));
+        }
+
+        // 生成随机最简化行阶梯矩阵 前rank_col中选择nrows个主元列
+        void randomize(size_t rank_col, uint_fast32_t seed)
+        {
+            assert(nrows <= rank_col && rank_col <= ncols);
+            randomize(seed);
+            vector<size_t> pivot(rank_col);
+            iota(pivot.begin(), pivot.end(), 0);
+            random_shuffle(pivot.begin(), pivot.end());
+            pivot.resize(nrows);
+            sort(pivot.begin(), pivot.end());
+
+            vector<bool> pivotmask(width, true);
+            for (size_t r = 0; r < nrows; r++)
+            {
+                pivotmask[pivot[r]] = false;
+            }
+
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < pivot[r]; w++)
+                {
+                    *at_base(r, w) = base_zero;
+                }
+                *at_base(r, pivot[r]) = base_one;
+                for (size_t w = pivot[r] + 1; w < rank_col; w++)
+                {
+                    if (!pivotmask[w])
+                    {
+                        *at_base(r, w) = base_zero;
+                    }
+                }
+            }
+        }
+
+        bool operator==(const MatGFP &m) const
+        {
+            if (nrows != m.nrows || ncols != m.ncols)
+            {
+                return false;
+            }
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    if (*at_base(r, w) != *m.at_base(r, w))
+                    {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+
+        bool operator==(const gfp_t base) const
+        {
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    if (*at_base(r, w) != base)
+                    {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+        void operator+=(const MatGFP &m)
+        {
+            assert(nrows == m.nrows && ncols == m.ncols);
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    *at_base(r, w) = (*at_base(r, w) + *m.at_base(r, w)) % gfprime;
+                }
+            }
+        }
+        MatGFP operator+(const MatGFP &m) const
+        {
+            MatGFP temp(*this);
+            temp += m;
+            return temp;
+        }
+        // a(m*k)*b(k,n) 其中k不超过65536
+        void cpu_addmul(const MatGFP &a, const MatGFP &b)
+        {
+            assert(a.ncols == b.nrows && a.nrows == nrows && b.ncols == ncols);
+            assert(a.ncols <= 65536);
+            for (size_t r = 0; r < nrows; r++)
+            {
+                for (size_t w = 0; w < width; w++)
+                {
+                    for (size_t i = 0; i < a.ncols; i++)
+                    {
+                        *at_base(r, w) += (*a.at_base(r, i) * *b.at_base(i, w)) % gfprime;
+                    }
+                    *at_base(r, w) %= gfprime;
+                }
+            }
+        }
+        void gpu_mul(const MatGFP &a, const MatGFP &b);
+
+        MatGFP operator*(const MatGFP &m) const
+        {
+            MatGFP temp(nrows, m.width);
+            temp.gpu_mul(*this, m);
+            return temp;
+        }
+
+        void cpu_swap_row(size_t r1, size_t r2)
+        {
+            if (r1 == r2)
+            {
+                return;
+            }
+            gfp_t *p1 = at_base(r1, 0);
+            gfp_t *p2 = at_base(r2, 0);
+            for (size_t i = 0; i < width; i++)
+            {
+                gfp_t temp = p1[i];
+                p1[i] = p2[i];
+                p2[i] = temp;
+            }
+        }
+
+        void cpu_mul_row(size_t r, gfp_t mul, size_t begin_w = 0, size_t end_w = 0)
+        {
+            for (size_t i = begin_w; i < (end_w ? end_w : width); i++)
+            {
+                *at_base(r, i) = (*at_base(r, i) * mul) % gfprime;
+            }
+        }
+
+        void cpu_addmul_row(size_t r1, size_t r2, gfp_t mul, size_t begin_w = 0, size_t end_w = 0)
+        {
+            for (size_t i = begin_w; i < (end_w ? end_w : width); i++)
+            {
+                *at_base(r1, i) = (*at_base(r1, i) + *at_base(r2, i) * mul) % gfprime;
+            }
+        }
+
+        ElimResult gpu_elim();
+        ElimResult cpu_elim();
+
+        friend ostream &operator<<(ostream &out, const MatGFP &m);
+
+        size_t nrows, ncols, width;
+
+    private:
+        MatGFP() : nrows(0), ncols(0), width(0), rowstride(0), type(moved), data(nullptr) {}
+
+        size_t rowstride;
+        MatType type;
+        gfp_t *data;
+    };
+
+    ostream &operator<<(ostream &out, const MatGFP &m)
+    {
+        for (size_t r = 0; r < m.nrows; r++)
+        {
+            for (size_t w = 0; w < m.width; w++)
+            {
+                printf("%05u ", *m.at_base(r, w));
+            }
+            printf("\n");
+        }
+        return out;
+    }
+}
+
+#endif

include/gfp/gfp_mul.cuh

@@ -0,0 +1,85 @@
+#ifndef GFP_MUL_CUH
+#define GFP_MUL_CUH
+
+#include "gfp_mat.cuh"
+
+namespace gfp
+{
+
+    static const int BlockRow = 128, BlockCol = 128; // 每个block处理c矩阵的一个子块
+    static const int StepSize = 8;                   // block中一个循环处理的A矩阵的列数（B矩阵的行数）
+
+    static_assert(BlockCol % THREAD_X == 0 && BlockRow % THREAD_Y == 0);
+
+    __global__ void gpu_mul_kernel(gfp_t *__restrict__ a, const size_t a_rs, gfp_t *__restrict__ b, const size_t b_rs, gfp_t *__restrict__ c, const size_t c_rs, const size_t nrows, const size_t ncols, const size_t nsteps)
+    {
+
+        const unsigned int bx = blockIdx.x;
+        const unsigned int by = blockIdx.y;
+        const unsigned int tx = threadIdx.x;
+        const unsigned int ty = threadIdx.y;
+        const unsigned int tid = ty * blockDim.x + tx;
+
+        __shared__ gfp_t s_a[StepSize][BlockRow];
+        __shared__ gfp_t s_b[StepSize][BlockCol];
+
+        gfp_t tmp_c[BlockRow / THREAD_Y][BlockCol / THREAD_X] = {0};
+
+        for (int s = 0; s < (nsteps - 1) / StepSize + 1; s++)
+        {
+            for (int k = tid; k < StepSize * BlockRow; k += blockDim.x * blockDim.y)
+            {
+                const int a_r = k / StepSize;
+                const int a_c = k % StepSize;
+                s_a[a_c][a_r] = by * BlockRow + a_r < nrows && s * StepSize + a_c < nsteps ? *at_base(a, a_rs, by * BlockRow + a_r, s * StepSize + a_c) : 0;
+                const int b_r = k / BlockCol;
+                const int b_c = k % BlockCol;
+                s_b[b_r][b_c] = s * StepSize + b_r < nsteps && bx * BlockCol + b_c < ncols ? *at_base(b, b_rs, s * StepSize + b_r, bx * BlockCol + b_c) : 0;
+            }
+            __syncthreads();
+            for (int k = 0; k < StepSize; k++)
+            {
+                for (int j = 0; j < BlockRow / THREAD_Y; j++)
+                {
+                    for (int i = 0; i < BlockCol / THREAD_X; i++)
+                    {
+                        tmp_c[j][i] += (s_a[k][j * THREAD_Y + ty] * s_b[k][i * THREAD_X + tx]) % gfprime;
+                    }
+                }
+            }
+            __syncthreads();
+            if (s & gfp_fullmask == gfp_fullmask)
+            {
+                for (int j = 0; j < BlockRow / THREAD_Y; j++)
+                {
+                    for (int i = 0; i < BlockCol / THREAD_X; i++)
+                    {
+                        tmp_c[j][i] %= gfprime;
+                    }
+                }
+            }
+        }
+        for (int j = 0; j < BlockRow / THREAD_Y; j++)
+        {
+            for (int i = 0; i < BlockCol / THREAD_X; i++)
+            {
+                if (by * BlockRow + j * THREAD_Y + ty < nrows && bx * BlockCol + i * THREAD_X + tx < ncols)
+                {
+                    *at_base(c, c_rs, by * BlockRow + j * THREAD_Y + ty, bx * BlockCol + i * THREAD_X + tx) = tmp_c[j][i] % gfprime;
+                }
+            }
+        }
+    }
+
+    __host__ void MatGFP::gpu_mul(const MatGFP &a, const MatGFP &b)
+    {
+        assert(a.ncols == b.nrows && a.nrows == nrows && b.ncols == ncols);
+
+        dim3 block(THREAD_X, THREAD_Y);
+        dim3 grid((width - 1) / block.x + 1, (nrows - 1) / block.y + 1);
+        gpu_mul_kernel<<<grid, block>>>(a.data, a.rowstride, b.data, b.rowstride, data, rowstride, nrows, width, a.width);
+        cudaDeviceSynchronize();
+    }
+}
+
+#endif

include/header.cuh

@@ -3,89 +3,27 @@
 
 #include <iostream>
 #include <cassert>
-#include <cuda_runtime.h>
 
-#include <set> // gf28
-// #include <map>
-// #include <vector>
-
-// #include <algorithm>
-// #include <numeric>
-#include <random> // matrix
-// #include <omp.h>
+#include <cpp_progress.hpp>
 
 using namespace std;
 
 using base_t = uint64_t;
-using gf28_t = uint8_t;
 
-static const size_t base_deg = 8;
-static const size_t base_num = 8;
-static const size_t base_len = 64;
-static_assert(base_len == base_deg * base_num && base_len == sizeof(base_t) * 8);
+static const size_t base_len = sizeof(base_t) * 8;
 
 static const base_t base_zero = (base_t)0x00'00'00'00'00'00'00'00;
 static const base_t base_one = (base_t)0x00'00'00'00'00'00'00'01;
-static const gf28_t gf28_zero = (gf28_t)0x00;
-static const gf28_t gf28_one = (gf28_t)0x01;
 
 static const base_t base_fullmask = (base_t)0xFF'FF'FF'FF'FF'FF'FF'FF;
-static const base_t base_deg_mask[8] = {
-    (base_t)0x00'00'00'00'00'00'00'FF,
-    (base_t)0x00'00'00'00'00'00'FF'00,
-    (base_t)0x00'00'00'00'00'FF'00'00,
-    (base_t)0x00'00'00'00'FF'00'00'00,
-    (base_t)0x00'00'00'FF'00'00'00'00,
-    (base_t)0x00'00'FF'00'00'00'00'00,
-    (base_t)0x00'FF'00'00'00'00'00'00,
-    (base_t)0xFF'00'00'00'00'00'00'00};
 
-static const size_t THREAD_X = 32;       // 列
-static const size_t THREAD_Y = base_deg; // 行
+static const size_t THREAD_X = 16; // 列
+static const size_t THREAD_Y = 16; // 行
 
-__constant__ gf28_t d_mul_table[1 << base_deg][1 << base_deg];
-
-__device__ base_t *at_pitch(base_t *base, size_t pitch, size_t r, size_t w)
+template <typename T>
+__host__ __device__ T *at_base(T *base, size_t rowstride, size_t r, size_t w)
 {
-    return base + r * pitch + w;
-}
-
-__host__ __device__ inline size_t offset(size_t idx)
-{
-    return idx << 3;
-}
-
-__host__ __device__ inline gf28_t get8(base_t src, size_t idx)
-{
-    return (gf28_t)(src >> offset(idx));
-}
-
-// 确保set8对应位置的值为0
-__host__ __device__ inline void set8(base_t &dst, size_t idx, gf28_t src)
-{
-    dst |= (base_t)src << offset(idx);
-}
-
-__host__ inline void del8(base_t &dst, size_t idx)
-{
-    dst &= ~base_deg_mask[idx];
-}
-
-__device__ base_t mul_base(const gf28_t val, const base_t base, const size_t offset = 0)
-{
-    base_t temp = base_zero;
-    for (size_t i = offset; i < base_num; i++)
-    {
-        set8(temp, i, d_mul_table[val][get8(base, i)]);
-    }
-    return temp;
-}
-
-__host__ inline base_t rev8(base_t n)
-{
-    n = (n & 0xff00ff00ff00ff00ul) >> 8 | (n & 0x00ff00ff00ff00fful) << 8;
-    n = (n & 0xffff0000ffff0000ul) >> 16 | (n & 0x0000ffff0000fffful) << 16;
-    return n >> 32 | n << 32;
+    return base + r * rowstride + w;
 }
 
 #define CUDA_CHECK(call)                                                 \

include/matrix.cuh

@@ -1,192 +0,0 @@
-#ifndef MATRIX_CUH
-#define MATRIX_CUH
-
-#include "header.cuh"
-#include "gf28.cuh"
-
-class GF28Matrix
-{
-public:
-    enum MatType
-    {
-        root,
-        view
-    };
-    // 只能构造root矩阵
-    GF28Matrix(size_t nrows, size_t ncols) : nrows(nrows), ncols(ncols), type(root)
-    {
-        width = (ncols - 1) / base_num + 1;
-        pitch = ((width - 1) / 4) * 4 + 1; // 以32字节（4*64bit）为单位对齐
-        CUDA_CHECK(cudaMallocManaged((void **)&data, nrows * pitch * sizeof(base_t)));
-        CUDA_CHECK(cudaMemset(data, 0, nrows * pitch * sizeof(base_t)));
-    }
-    // 只能拷贝构造root矩阵
-    GF28Matrix(const GF28Matrix &m) : GF28Matrix(m.nrows, m.ncols)
-    {
-        cudaMemcpy2D(data, pitch * sizeof(base_t), m.data, m.pitch * sizeof(base_t), m.width * sizeof(base_t), nrows, cudaMemcpyDefault);
-    }
-    GF28Matrix(GF28Matrix &&m) noexcept : nrows(m.nrows), ncols(m.ncols), width(m.width), pitch(m.pitch), type(m.type), data(m.data)
-    {
-        m.nrows = 0;
-        m.ncols = 0;
-        m.width = 0;
-        m.pitch = 0;
-        m.type = view;
-        m.data = nullptr;
-    }
-    GF28Matrix &operator=(const GF28Matrix &m)
-    {
-        if (this == &m)
-        {
-            return *this;
-        }
-        assert(nrows == m.nrows && ncols == m.ncols);
-        cudaMemcpy2D(data, pitch * sizeof(base_t), m.data, m.pitch * sizeof(base_t), m.width * sizeof(base_t), nrows, cudaMemcpyDefault);
-        return *this;
-    }
-    GF28Matrix &operator=(GF28Matrix &&m) noexcept
-    {
-        if (this == &m)
-        {
-            return *this;
-        }
-        if (type == root)
-        {
-            CUDA_CHECK(cudaFree(data));
-        }
-        nrows = m.nrows;
-        ncols = m.ncols;
-        width = m.width;
-        pitch = m.pitch;
-        type = m.type;
-        data = m.data;
-        m.nrows = 0;
-        m.ncols = 0;
-        m.width = 0;
-        m.pitch = 0;
-        m.type = view;
-        m.data = nullptr;
-        return *this;
-    }
-
-    ~GF28Matrix()
-    {
-        if (type == root)
-        {
-            CUDA_CHECK(cudaFree(data));
-        }
-    }
-
-    inline base_t *at_base(size_t r, size_t w) const
-    {
-        return data + r * pitch + w;
-    }
-
-    // 只能以base_t为单位进行操作
-    GF28Matrix createView(size_t begin_ri, size_t begin_wi, size_t end_rj, size_t end_wj) const
-    {
-        assert(begin_ri < end_rj && end_rj <= nrows && begin_wi < end_wj && end_wj <= width);
-        GF28Matrix view;
-        view.nrows = end_rj - begin_ri;
-        view.ncols = (end_wj == width ? ncols : end_wj * base_num) - begin_wi * base_num;
-        view.width = end_wj - begin_wi;
-        view.pitch = pitch;
-        view.data = at_base(begin_ri, begin_wi);
-        return view;
-    }
-
-    void randomize(base_t seed)
-    {
-        assert(type == root);
-        static default_random_engine e(seed);
-        static uniform_int_distribution<base_t> d;
-        base_t lastmask = base_fullmask >> (width * base_len - ncols * base_deg);
-        for (size_t r = 0; r < nrows; r++)
-        {
-            for (size_t w = 0; w < width; w++)
-            {
-                *at_base(r, w) = d(e);
-            }
-            *at_base(r, width - 1) &= lastmask;
-        }
-    }
-
-    bool operator==(const GF28Matrix &m) const
-    {
-        if (nrows != m.nrows || ncols != m.ncols)
-        {
-            return false;
-        }
-        for (size_t r = 0; r < nrows; r++)
-        {
-            for (size_t w = 0; w < width; w++)
-            {
-                if (*at_base(r, w) != *m.at_base(r, w))
-                {
-                    return false;
-                }
-            }
-        }
-        return true;
-    }
-
-    bool operator==(const base_t base) const
-    {
-        for (size_t r = 0; r < nrows; r++)
-        {
-            for (size_t w = 0; w < width; w++)
-            {
-                if (*at_base(r, w) != base)
-                {
-                    return false;
-                }
-            }
-        }
-        return true;
-    }
-    void operator^=(const GF28Matrix &m)
-    {
-        assert(nrows == m.nrows && ncols == m.ncols);
-        for (size_t r = 0; r < nrows; r++)
-        {
-            for (size_t w = 0; w < width; w++)
-            {
-                *at_base(r, w) ^= *m.at_base(r, w);
-            }
-        }
-    }
-    GF28Matrix operator^(const GF28Matrix &m) const
-    {
-        GF28Matrix temp(*this);
-        temp ^= m;
-        return temp;
-    }
-
-    friend ostream &operator<<(ostream &out, const GF28Matrix &m);
-    void gpu_addmul(const GF28Matrix &a, const GF28Matrix &b, const GF28 &gf);
-
-    // size_t nrows, ncols;
-    // size_t width, pitch;
-
-private:
-    GF28Matrix() : nrows(0), ncols(0), width(0), pitch(0), type(view), data(nullptr) {}
-    size_t nrows, ncols;
-    size_t width, pitch;
-    MatType type;
-    base_t *data;
-};
-
-ostream &operator<<(ostream &out, const GF28Matrix &m)
-{
-    for (size_t r = 0; r < m.nrows; r++)
-    {
-        for (size_t w = 0; w < m.width; w++)
-        {
-            printf("%016lX ", rev8(*m.at_base(r, w)));
-        }
-        printf("\n");
-    }
-    return out;
-}
-
-#endif

include/multiplication.cuh

@@ -1,52 +0,0 @@
-#ifndef MULTIPLICATION_CUH
-#define MULTIPLICATION_CUH
-
-#include "matrix.cuh"
-#include "gf28.cuh"
-
-// 处理32base列的所有行
-__global__ void gpu_addmul_kernel(base_t *a, size_t a_pitch, base_t *b, size_t b_pitch, base_t *c, size_t c_pitch, size_t nrows, size_t width)
-{
-    __shared__ __align__(8) base_t src[base_deg][THREAD_X];
-    size_t r = threadIdx.y;
-    size_t w = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (w >= width)
-    {
-        return;
-    }
-
-    if (r < nrows && w < width)
-        src[threadIdx.y][threadIdx.x] = *at_pitch(b, b_pitch, r, w);
-    else
-        src[threadIdx.y][threadIdx.x] = base_zero;
-
-    __syncthreads();
-
-    for (; r < nrows; r += base_deg)
-    {
-        base_t val = *at_pitch(a, a_pitch, r, 0);
-        base_t temp = base_zero;
-        for (size_t i = 0; i < base_deg; i++)
-        {
-            temp ^= mul_base(get8(val, i), src[i][threadIdx.x]);
-        }
-        *at_pitch(c, c_pitch, r, w) ^= temp;
-    }
-}
-
-__host__ void GF28Matrix::gpu_addmul(const GF28Matrix &a, const GF28Matrix &b, const GF28 &gf)
-{
-    assert(a.ncols == b.nrows && a.nrows == nrows && b.ncols == ncols);
-    cudaMemcpyToSymbol(d_mul_table, gf.mul_table, (1 << base_deg) * (1 << base_deg) * sizeof(gf28_t));
-    for (size_t w = 0; w < a.width; w++)
-    {
-        dim3 block(THREAD_X, THREAD_Y);
-        dim3 grid((b.width - 1) / block.x + 1);
-        cudaDeviceSynchronize();
-        gpu_addmul_kernel<<<grid, block>>>(a.at_base(0, w), a.pitch, b.at_base(w * base_num, 0), b.pitch, data, pitch, nrows, width);
-    }
-    cudaDeviceSynchronize();
-}
-
-#endif

src/main.cu

@@ -1,11 +1,30 @@
+// #define SHOW_PROGRESS_BAR
+
 #include "cuelim.cuh"
+// #include <m4ri/m4ri.h>
+
+// #include "interface.cuh"
+
+// #undef SHOW_PROGRESS_BAR
+
+using namespace gf2;
+
+bool test_gf2_elim(size_t rank, size_t rank_col, size_t nrows, size_t ncols, uint_fast32_t seed)
+{
+    assert(rank <= nrows && rank <= rank_col && rank_col <= ncols);
+    MatGF2 rdc(rank, ncols);
+    rdc.randomize(rank_col, seed);
+    MatGF2 mix(nrows, rank);
+    mix.randomize(seed);
+    MatGF2 src = mix * rdc;
+
+    ElimResult res = src.gpu_elim();
+
+    MatGF2 win(src, 0, 0, res.rank, src.width);
+    return rdc == win;
+}
 
 int main()
 {
-    GF28Matrix a(10000, 10000);
-    a.randomize(123);
-    GF28Matrix b(10000, 10000);
-    b.randomize(456);
-    GF28Matrix c(10000, 10000);
-    c.gpu_addmul(a, b, GF28(0b100011101));
+    cout << test_gf2_elim(480, 960, 600, 1200, 123) << endl;
 }

test/CMakeLists.txt

@@ -3,8 +3,12 @@ find_package(GTest REQUIRED)  # 查找GTest库
 include_directories(${PROJECT_SOURCE_DIR}/test)  # 添加测试头文件目录
 
 set(TEST_SRC_FILES
-    "test_gf28.cu"
-    "test_matrix.cu"
+    "test_gf256_header.cu"
+    "test_gf256_matrix.cu"
+    "test_gf256_elim.cu"
+    "test_gfp_mul.cu"
+    "test_gfp_elim.cu"
+    "test_gf2_elim.cu"
 )
 
 foreach(SRC ${TEST_SRC_FILES})
@@ -14,13 +18,14 @@ foreach(SRC ${TEST_SRC_FILES})
     gtest_discover_tests(${SRC_NAME})
 endforeach()
 
-# set(TEST_M4RIE_SRC_FILES
-#     "test_m4rie_interface.cu"
-# )
+set(TEST_M4RIE_SRC_FILES
+    "test_m4ri_interface.cu"
+    "test_m4rie_interface.cu"
+)
 
-# foreach(SRC ${TEST_M4RIE_SRC_FILES})
-#     get_filename_component(SRC_NAME ${SRC} NAME_WE)
-#     add_executable(${SRC_NAME} ${SRC})
-#     target_link_libraries(${SRC_NAME} GTest::GTest GTest::Main m4ri m4rie)
-#     gtest_discover_tests(${SRC_NAME})
-# endforeach()
+foreach(SRC ${TEST_M4RIE_SRC_FILES})
+    get_filename_component(SRC_NAME ${SRC} NAME_WE)
+    add_executable(${SRC_NAME} ${SRC})
+    target_link_libraries(${SRC_NAME} GTest::GTest GTest::Main m4ri m4rie)
+    gtest_discover_tests(${SRC_NAME})
+endforeach()

test/test_gf256_elim.cu

@@ -0,0 +1,32 @@
+#include <gtest/gtest.h>
+#include "test_header.cuh"
+
+using namespace gf256;
+
+bool test_gf256_elim(size_t rank, size_t rank_col, size_t nrows, size_t ncols, const GF256 &gf256, uint_fast32_t seed)
+{
+    assert(rank <= nrows && rank <= rank_col && rank_col <= ncols);
+    MatGF256 rdc(rank, ncols);
+    rdc.randomize(rank_col, seed);
+    MatGF256 mix(nrows, rank);
+    mix.randomize(seed);
+    MatGF256 src = gpu_mul(mix, rdc, gf256);
+    ElimResult res = src.gpu_elim(gf256);
+    MatGF256 win(src, 0, 0, res.rank, src.width);
+    return rdc == win;
+}
+
+TEST(TestGF256Elim, Small)
+{
+    uint_fast32_t seed = 41921095;
+    GF256 gf256(0b100011101);
+    EXPECT_TRUE(test_gf256_elim(5, 7, 6, 8, gf256, seed));
+}
+
+TEST(TestGF256Elim, Mediem)
+{
+    uint_fast32_t seed = 41921095;
+    GF256 gf256(0b100011101);
+    EXPECT_TRUE(test_gf256_elim(50, 70, 60, 80, gf256, seed));
+    EXPECT_TRUE(test_gf256_elim(500, 700, 600, 800, gf256, seed));
+}

test/test_gf256_header.cu

@@ -1,7 +1,9 @@
 #include <gtest/gtest.h>
 #include "test_header.cuh"
 
-vector<gf28_t> expect_inv_table{
+using namespace gf256;
+
+vector<gf256_t> expect_inv_table{
     0x00, 0x01, 0x8E, 0xF4, 0x47, 0xA7, 0x7A, 0xBA, 0xAD, 0x9D, 0xDD, 0x98, 0x3D, 0xAA, 0x5D, 0x96,
     0xD8, 0x72, 0xC0, 0x58, 0xE0, 0x3E, 0x4C, 0x66, 0x90, 0xDE, 0x55, 0x80, 0xA0, 0x83, 0x4B, 0x2A,
     0x6C, 0xED, 0x39, 0x51, 0x60, 0x56, 0x2C, 0x8A, 0x70, 0xD0, 0x1F, 0x4A, 0x26, 0x8B, 0x33, 0x6E,
@@ -19,11 +21,11 @@ vector<gf28_t> expect_inv_table{
     0x14, 0x3F, 0xE6, 0xF0, 0x86, 0xB1, 0xE2, 0xF1, 0xFA, 0x74, 0xF3, 0xB4, 0x6D, 0x21, 0xB2, 0x6A,
     0xE3, 0xE7, 0xB5, 0xEA, 0x03, 0x8F, 0xD3, 0xC9, 0x42, 0xD4, 0xE8, 0x75, 0x7F, 0xFF, 0x7E, 0xFD};
 
-TEST(TestGF28, Inv)
+TEST(TestGF256Header, Inv)
 {
-    GF28 gf28(0b100011101);
-    for (size_t x = 0; x < 1 << base_deg; x++)
+    GF256 gf256(0b100011101);
+    for (size_t x = 0; x < 1 << gf256_len; x++)
     {
-        EXPECT_EQ(gf28.inv(x), expect_inv_table[x]);
+        EXPECT_EQ(gf256.inv(x), expect_inv_table[x]);
     }
 }

test/test_gf256_matrix.cu

@@ -0,0 +1,33 @@
+#include <gtest/gtest.h>
+#include "test_header.cuh"
+
+using namespace gf256;
+
+TEST(TestGF256Matrix, Equal)
+{
+    MatGF256 a(50, 50);
+    EXPECT_TRUE(a == base_zero);
+    MatGF256 v(a, 0, 0, 30, 3);
+    EXPECT_TRUE(v == base_zero);
+    a.randomize(1243);
+    EXPECT_TRUE(a == a);
+    EXPECT_TRUE(v == v);
+    MatGF256 b(50, 50);
+    b.randomize(1243);
+    EXPECT_FALSE(a == b);
+}
+
+TEST(TestGF256Matrix, Xor)
+{
+    MatGF256 a(50, 50), b(50, 50);
+    a.randomize(1243);
+    b.randomize(1243);
+    MatGF256 c = a ^ b;
+    a ^= c;
+    EXPECT_TRUE(a == b);
+    MatGF256 va(a, 20, 1, 30, 3);
+    MatGF256 vb(b, 10, 2, 20, 4);
+    MatGF256 vc = va ^ vb;
+    va ^= vc;
+    EXPECT_TRUE(va == vb);
+}

test/test_gf2_elim.cu

@@ -0,0 +1,30 @@
+#include <gtest/gtest.h>
+#include "test_header.cuh"
+
+using namespace gf2;
+
+bool test_gf2_elim(size_t rank, size_t rank_col, size_t nrows, size_t ncols, uint_fast32_t seed)
+{
+    assert(rank <= nrows && rank <= rank_col && rank_col <= ncols);
+    MatGF2 rdc(rank, ncols);
+    rdc.randomize(rank_col, seed);
+    MatGF2 mix(nrows, rank);
+    mix.randomize(seed);
+    MatGF2 src = mix * rdc;
+    ElimResult res = src.gpu_elim();
+    MatGF2 win(src, 0, 0, res.rank, src.width);
+    return rdc == win;
+}
+
+TEST(TestGF2Elim, Small)
+{
+    uint_fast32_t seed = 41921095;
+    EXPECT_TRUE(test_gf2_elim(5, 7, 6, 8, seed));
+}
+
+TEST(TestGF2Elim, Mediem)
+{
+    uint_fast32_t seed = 41921095;
+    EXPECT_TRUE(test_gf2_elim(50, 70, 60, 80, seed));
+    EXPECT_TRUE(test_gf2_elim(500, 700, 600, 800, seed));
+}

test/test_gfp_elim.cu

@@ -0,0 +1,29 @@
+#include <gtest/gtest.h>
+#include "test_header.cuh"
+
+using namespace gfp;
+
+bool test_gfp_elim(size_t rank, size_t rank_col, size_t nrows, size_t ncols, uint_fast32_t seed)
+{
+    MatGFP rdc(rank, ncols);
+    rdc.randomize(rank_col, seed);
+    MatGFP mix(nrows, rank);
+    mix.randomize(seed);
+    MatGFP a = mix * rdc;
+    ElimResult res = a.gpu_elim();
+    MatGFP win(a, 0, 0, res.rank, a.width);
+    return rdc == win;
+}
+
+TEST(TestGFPMul, Small)
+{
+    uint_fast32_t seed = 41921095;
+    EXPECT_TRUE(test_gfp_elim(5, 6, 7, 8, seed));
+}
+
+TEST(TestGFPMul, Mediem)
+{
+    uint_fast32_t seed = 41921095;
+    EXPECT_TRUE(test_gfp_elim(50, 60, 70, 80, seed));
+    EXPECT_TRUE(test_gfp_elim(500, 600, 700, 800, seed));
+}

test/test_gfp_mul.cu

@@ -0,0 +1,29 @@
+#include <gtest/gtest.h>
+#include "test_header.cuh"
+
+using namespace gfp;
+
+bool test_gfp_mul(size_t m, size_t k, size_t n, uint_fast32_t seed)
+{
+    MatGFP a(m, k);
+    MatGFP b(k, n);
+    MatGFP c(m, n);
+    a.randomize(seed);
+    b.randomize(seed);
+    c.cpu_addmul(a, b);
+    MatGFP d = a * b;
+    return c == d;
+}
+
+TEST(TestGFPMul, Small)
+{
+    uint_fast32_t seed = 41921095;
+    EXPECT_TRUE(test_gfp_mul(5, 7, 6, seed));
+}
+
+TEST(TestGFPMul, Mediem)
+{
+    uint_fast32_t seed = 41921095;
+    EXPECT_TRUE(test_gfp_mul(50, 70, 60, seed));
+    EXPECT_TRUE(test_gfp_mul(500, 700, 600, seed));
+}

test/test_m4ri_interface.cu

@@ -0,0 +1,37 @@
+#include <gtest/gtest.h>
+#include "test_header.cuh"
+#include "cuelim_m4ri.cuh"
+
+using namespace gf2;
+
+bool test_gf2_elim_interface(size_t rank, size_t rank_col, size_t nrows, size_t ncols, uint_fast32_t seed)
+{
+    assert(rank <= nrows && rank <= rank_col && rank_col <= ncols);
+    MatGF2 rdc(rank, ncols);
+    rdc.randomize(rank_col, seed);
+    MatGF2 mix(nrows, rank);
+    mix.randomize(seed);
+    MatGF2 src = mix * rdc;
+
+    mzd_t *A_m4ri = mzd_init(src.nrows, src.ncols);
+    mzdwrite(src, A_m4ri);
+    mzd_t *A_m4ri_copy = mzd_copy(NULL, A_m4ri);
+
+    base_t rank_interface = gpu_mzd_elim(A_m4ri);
+    rci_t rank_m4rie = mzd_echelonize_m4ri(A_m4ri_copy, 1, 8);
+
+    return (rank_interface == rank_m4rie) && (mzd_cmp(A_m4ri, A_m4ri_copy) == 0);
+}
+
+TEST(TestM4riInterface, Small)
+{
+    uint_fast32_t seed = 41921095;
+    EXPECT_TRUE(test_gf2_elim_interface(5, 7, 6, 8, seed));
+}
+
+TEST(TestM4riInterface, Mediem)
+{
+    uint_fast32_t seed = 41921095;
+    EXPECT_TRUE(test_gf2_elim_interface(50, 70, 60, 80, seed));
+    EXPECT_TRUE(test_gf2_elim_interface(500, 700, 600, 800, seed));
+}

test/test_m4rie_interface.cu

@@ -0,0 +1,40 @@
+#include <gtest/gtest.h>
+#include "test_header.cuh"
+#include "cuelim_m4rie.cuh"
+
+using namespace gf256;
+
+bool test_gf256_elim_interface(size_t rank, size_t rank_col, size_t nrows, size_t ncols, const GF256 &gf256, uint_fast32_t seed)
+{
+    assert(rank <= nrows && rank <= rank_col && rank_col <= ncols);
+    MatGF256 rdc(rank, ncols);
+    rdc.randomize(rank_col, seed);
+    MatGF256 mix(nrows, rank);
+    mix.randomize(seed);
+    MatGF256 src = gpu_mul(mix, rdc, gf256);
+
+    gf2e *ff_m4rie = gf2e_init(gf256.poly());
+    mzed_t *A_m4rie = mzed_init(ff_m4rie, src.nrows, src.ncols);
+    mzedwrite(src, A_m4rie);
+    mzed_t *A_m4rie_copy = mzed_copy(NULL, A_m4rie);
+
+    base_t rank_interface = gpu_mzed_elim(A_m4rie);
+    rci_t rank_m4rie = mzed_echelonize_newton_john(A_m4rie_copy, 1);
+
+    return (rank_interface == rank_m4rie) && (mzed_cmp(A_m4rie, A_m4rie_copy) == 0);
+}
+
+TEST(TestM4rieInterface, Small)
+{
+    uint_fast32_t seed = 41921095;
+    GF256 gf256(0b100011101);
+    EXPECT_TRUE(test_gf256_elim_interface(5, 7, 6, 8, gf256, seed));
+}
+
+TEST(TestM4rieInterface, Mediem)
+{
+    uint_fast32_t seed = 41921095;
+    GF256 gf256(0b100011101);
+    EXPECT_TRUE(test_gf256_elim_interface(50, 70, 60, 80, gf256, seed));
+    EXPECT_TRUE(test_gf256_elim_interface(500, 700, 600, 800, gf256, seed));
+}

test/test_matrix.cu

@@ -1,39 +0,0 @@
-#include <gtest/gtest.h>
-#include "test_header.cuh"
-
-TEST(TestMatrix, Equal)
-{
-    GF28Matrix a(50, 50);
-    EXPECT_TRUE(a == base_zero);
-    GF28Matrix v = a.createView(0, 0, 30, 3);
-    EXPECT_TRUE(v == base_zero);
-    a.randomize(1243);
-    EXPECT_TRUE(a == a);
-    EXPECT_TRUE(v == v);
-    GF28Matrix b(50, 50);
-    b.randomize(1243);
-    EXPECT_FALSE(a == b);
-}
-
-TEST(TestMatrix, Xor)
-{
-    GF28Matrix a(50, 50), b(50, 50);
-    a.randomize(1243);
-    b.randomize(1243);
-    GF28Matrix c = a ^ b;
-    a ^= c;
-    EXPECT_TRUE(a == b);
-    GF28Matrix va = a.createView(20, 1, 30, 3);
-    GF28Matrix vb = b.createView(10, 2, 20, 4);
-    GF28Matrix vc = va ^ vb;
-    va ^= vc;
-    EXPECT_TRUE(va == vb);
-}
-
-// TEST(TestMatrix, Basic)
-// {
-//     GF28Matrix a(50, 50);
-//     GF28Matrix v = a.createView(0, 0, 30, 3);
-
-//     EXPECT_EQ(v.type, GF28Matrix::view);
-// }