// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// BoolMatrix.hpp
//
// ブール行列 ({0,1} 上の行列)
//
// 主な API:
//   BoolMatrix(rows, cols)       — ゼロ行列
//   BoolMatrix::identity(n)      — 単位行列
//   operator*, +, &              — 行列積 (OR-AND), 和 (OR), AND
//   transitiveClosure()          — 推移閉包 (Warshall)
//   reachability(from)           — 到達可能頂点集合
//   isReflexive, isSymmetric, isTransitive — 関係の性質

#ifndef CALX_LINALG_BOOL_MATRIX_HPP
#define CALX_LINALG_BOOL_MATRIX_HPP

#include <vector>
#include <cstdint>
#include <stdexcept>
#include <algorithm>

namespace calx {

/**
 * @brief ブール行列 (m × n, {0,1} 成分)
 */
class BoolMatrix {
public:
    BoolMatrix(std::size_t rows, std::size_t cols)
        : rows_(rows), cols_(cols), data_(rows, std::vector<bool>(cols, false)) {}

    /// ビットベクトルから構築
    BoolMatrix(std::size_t rows, std::size_t cols,
               const std::vector<std::vector<bool>>& data)
        : rows_(rows), cols_(cols), data_(data) {
        data_.resize(rows);
        for (auto& row : data_) row.resize(cols, false);
    }

    static BoolMatrix identity(std::size_t n) {
        BoolMatrix m(n, n);
        for (std::size_t i = 0; i < n; ++i) m.data_[i][i] = true;
        return m;
    }

    std::size_t rows() const { return rows_; }
    std::size_t cols() const { return cols_; }

    bool get(std::size_t i, std::size_t j) const { return data_[i][j]; }
    void set(std::size_t i, std::size_t j, bool v) { data_[i][j] = v; }

    /// ブール行列積 (OR-AND): C[i][j] = OR_k (A[i][k] AND B[k][j])
    friend BoolMatrix operator*(const BoolMatrix& a, const BoolMatrix& b) {
        if (a.cols_ != b.rows_)
            throw std::invalid_argument("BoolMatrix: dimension mismatch");
        BoolMatrix c(a.rows_, b.cols_);
        for (std::size_t i = 0; i < a.rows_; ++i)
            for (std::size_t k = 0; k < a.cols_; ++k)
                if (a.data_[i][k])
                    for (std::size_t j = 0; j < b.cols_; ++j)
                        if (b.data_[k][j])
                            c.data_[i][j] = true;
        return c;
    }

    /// OR (成分ごと)
    friend BoolMatrix operator|(const BoolMatrix& a, const BoolMatrix& b) {
        checkSameSize(a, b);
        BoolMatrix c(a.rows_, a.cols_);
        for (std::size_t i = 0; i < a.rows_; ++i)
            for (std::size_t j = 0; j < a.cols_; ++j)
                c.data_[i][j] = a.data_[i][j] || b.data_[i][j];
        return c;
    }

    /// AND (成分ごと)
    friend BoolMatrix operator&(const BoolMatrix& a, const BoolMatrix& b) {
        checkSameSize(a, b);
        BoolMatrix c(a.rows_, a.cols_);
        for (std::size_t i = 0; i < a.rows_; ++i)
            for (std::size_t j = 0; j < a.cols_; ++j)
                c.data_[i][j] = a.data_[i][j] && b.data_[i][j];
        return c;
    }

    /// 転置
    BoolMatrix transpose() const {
        BoolMatrix t(cols_, rows_);
        for (std::size_t i = 0; i < rows_; ++i)
            for (std::size_t j = 0; j < cols_; ++j)
                t.data_[j][i] = data_[i][j];
        return t;
    }

    /// 推移閉包 (Warshall のアルゴリズム)
    BoolMatrix transitiveClosure() const {
        if (rows_ != cols_)
            throw std::invalid_argument("BoolMatrix: transitive closure requires square matrix");
        BoolMatrix t = *this;
        for (std::size_t k = 0; k < rows_; ++k)
            for (std::size_t i = 0; i < rows_; ++i)
                if (t.data_[i][k])
                    for (std::size_t j = 0; j < rows_; ++j)
                        if (t.data_[k][j])
                            t.data_[i][j] = true;
        return t;
    }

    /// 頂点 from から到達可能な頂点集合
    std::vector<bool> reachability(std::size_t from) const {
        if (rows_ != cols_)
            throw std::invalid_argument("BoolMatrix: reachability requires square matrix");
        std::vector<bool> visited(rows_, false);
        std::vector<std::size_t> stack = {from};
        visited[from] = true;
        while (!stack.empty()) {
            auto v = stack.back(); stack.pop_back();
            for (std::size_t j = 0; j < cols_; ++j) {
                if (data_[v][j] && !visited[j]) {
                    visited[j] = true;
                    stack.push_back(j);
                }
            }
        }
        return visited;
    }

    /// 反射的か (対角成分がすべて 1)
    bool isReflexive() const {
        if (rows_ != cols_) return false;
        for (std::size_t i = 0; i < rows_; ++i)
            if (!data_[i][i]) return false;
        return true;
    }

    /// 対称か
    bool isSymmetric() const {
        if (rows_ != cols_) return false;
        for (std::size_t i = 0; i < rows_; ++i)
            for (std::size_t j = i + 1; j < cols_; ++j)
                if (data_[i][j] != data_[j][i]) return false;
        return true;
    }

    /// 推移的か (A² ⊆ A)
    bool isTransitive() const {
        if (rows_ != cols_) return false;
        auto sq = *this * *this;
        for (std::size_t i = 0; i < rows_; ++i)
            for (std::size_t j = 0; j < cols_; ++j)
                if (sq.data_[i][j] && !data_[i][j]) return false;
        return true;
    }

    friend bool operator==(const BoolMatrix& a, const BoolMatrix& b) {
        if (a.rows_ != b.rows_ || a.cols_ != b.cols_) return false;
        return a.data_ == b.data_;
    }

    friend bool operator!=(const BoolMatrix& a, const BoolMatrix& b) {
        return !(a == b);
    }

private:
    std::size_t rows_, cols_;
    std::vector<std::vector<bool>> data_;

    static void checkSameSize(const BoolMatrix& a, const BoolMatrix& b) {
        if (a.rows_ != b.rows_ || a.cols_ != b.cols_)
            throw std::invalid_argument("BoolMatrix: size mismatch");
    }
};

} // namespace calx

#endif // CALX_LINALG_BOOL_MATRIX_HPP