Skip to content

xor_basis.hpp

SECTIONMath INCLUDEnoya/xor_basis.hpp

Reduced linear basis over GF(2) for an unsigned integer type.

Verified by intersection_of_f2_vector_spaces.

\[ \displaystyle \operatorname{span}_{\mathbb F_2}(B)=\left\{\bigoplus_{b\in S}b:S\subseteq B\right\} \]

Implementation

View on GitHub

#ifndef NOYA_XOR_BASIS_HPP
#define NOYA_XOR_BASIS_HPP 1

/// @complexity Time: O(B^2) reduced insertion and O(B) membership/min/max query.
/// Space: O(B).

#include <array>
#include <cassert>
#include <cstdint>
#include <limits>
#include <type_traits>
#include <vector>

namespace noya {

/// @brief Reduced linear basis over GF(2) for an unsigned integer type.
template <class UInt, int bits = std::numeric_limits<UInt>::digits>
struct xor_basis {
  static_assert(std::is_integral_v<UInt> && std::is_unsigned_v<UInt>);
  static_assert(0 < bits && bits <= std::numeric_limits<UInt>::digits);

  std::array<UInt, bits> basis{};
  int dimension = 0;

  /// @brief Insert a value; return false if it is already in the span.
  bool insert(UInt value) {
    for (int bit = bits - 1; bit >= 0; bit--) {
      if (((value >> bit) & UInt(1)) == 0) {
        continue;
      }
      if (basis[bit] != 0) {
        value ^= basis[bit];
        continue;
      }
      basis[bit] = value;
      for (int lower = 0; lower < bit; lower++) {
        if ((basis[bit] >> lower) & UInt(1)) {
          basis[bit] ^= basis[lower];
        }
      }
      for (int higher = bit + 1; higher < bits; higher++) {
        if ((basis[higher] >> bit) & UInt(1)) {
          basis[higher] ^= basis[bit];
        }
      }
      dimension++;
      return true;
    }
    return false;
  }

  /// @brief Return whether value belongs to the represented xor span.
  bool contains(UInt value) const {
    for (int bit = bits - 1; bit >= 0; bit--) {
      if ((value >> bit) & UInt(1)) {
        value ^= basis[bit];
      }
    }
    return value == 0;
  }

  /// @brief Maximize seed xor x over all represented values x.
  UInt max_xor(UInt seed = 0) const {
    for (int bit = bits - 1; bit >= 0; bit--) {
      if ((seed ^ basis[bit]) > seed) {
        seed ^= basis[bit];
      }
    }
    return seed;
  }

  /// @brief Minimize seed xor x over all represented values x.
  UInt min_xor(UInt seed) const {
    for (int bit = bits - 1; bit >= 0; bit--) {
      if ((seed ^ basis[bit]) < seed) {
        seed ^= basis[bit];
      }
    }
    return seed;
  }

  /// @brief Return the k-th smallest distinct value in the span, including
  /// zero.
  UInt kth_smallest(std::uint64_t k) const {
    if (dimension < 64) {
      assert(k < (std::uint64_t(1) << dimension));
    }
    UInt result = 0;
    int index = 0;
    for (int bit = 0; bit < bits; bit++) {
      if (basis[bit] != 0) {
        if ((k >> index) & 1) {
          result ^= basis[bit];
        }
        index++;
      }
    }
    return result;
  }

  /// @brief Insert every basis vector from another span.
  void merge(const xor_basis &other) {
    for (UInt value : other.basis) {
      if (value != 0) {
        insert(value);
      }
    }
  }
};

/// @brief Return a basis of the intersection of two xor spans.  A homogeneous
/// system is built for U*a = V*b, with one equation per value bit. Gaussian
/// elimination produces a nullspace basis; projecting its U-coordinates gives
/// an independent basis of the intersection because both input lists are
/// required to be independent.
template <class UInt, int bits = std::numeric_limits<UInt>::digits>
std::vector<UInt>
xor_space_intersection(const std::vector<UInt> &first,
                       const std::vector<UInt> &second) {
  static_assert(std::is_integral_v<UInt> && std::is_unsigned_v<UInt>);
  static_assert(0 < bits && bits <= std::numeric_limits<UInt>::digits);
  const int left_size = int(first.size());
  const int variable_count = left_size + int(second.size());
  assert(variable_count <= 64);
  xor_basis<UInt, bits> left_check;
  xor_basis<UInt, bits> right_check;
  for (UInt value : first) {
    assert(left_check.insert(value));
  }
  for (UInt value : second) {
    assert(right_check.insert(value));
  }

  std::array<std::uint64_t, bits> equations{};
  for (int variable = 0; variable < variable_count; variable++) {
    UInt value = variable < left_size ? first[variable]
                                      : second[variable - left_size];
    for (int bit = 0; bit < bits; bit++) {
      if ((value >> bit) & UInt(1)) {
        equations[bit] |= std::uint64_t(1) << variable;
      }
    }
  }

  std::array<int, 64> pivot_row{};
  pivot_row.fill(-1);
  int rank = 0;
  for (int column = 0; column < variable_count; column++) {
    int pivot = rank;
    while (pivot < bits && ((equations[pivot] >> column) & 1) == 0) {
      pivot++;
    }
    if (pivot == bits) {
      continue;
    }
    std::swap(equations[rank], equations[pivot]);
    pivot_row[column] = rank;
    for (int row = 0; row < bits; row++) {
      if (row != rank && ((equations[row] >> column) & 1)) {
        equations[row] ^= equations[rank];
      }
    }
    rank++;
  }

  std::vector<UInt> result;
  for (int free_column = 0; free_column < variable_count; free_column++) {
    if (pivot_row[free_column] != -1) {
      continue;
    }
    std::uint64_t coefficients = std::uint64_t(1) << free_column;
    for (int column = 0; column < variable_count; column++) {
      int row = pivot_row[column];
      if (row != -1 && ((equations[row] >> free_column) & 1)) {
        coefficients |= std::uint64_t(1) << column;
      }
    }
    UInt value = 0;
    for (int index = 0; index < left_size; index++) {
      if ((coefficients >> index) & 1) {
        value ^= first[index];
      }
    }
    result.push_back(value);
  }
  return result;
}

} // namespace noya

#endif // NOYA_XOR_BASIS_HPP
#include <array>
#include <cassert>
#include <cstdint>
#include <limits>
#include <type_traits>
#include <vector>

/// @complexity Time: O(B^2) reduced insertion and O(B) membership/min/max query.
/// Space: O(B).

namespace noya {

/// @brief Reduced linear basis over GF(2) for an unsigned integer type.
template <class UInt, int bits = std::numeric_limits<UInt>::digits>
struct xor_basis {
  static_assert(std::is_integral_v<UInt> && std::is_unsigned_v<UInt>);
  static_assert(0 < bits && bits <= std::numeric_limits<UInt>::digits);

  std::array<UInt, bits> basis{};
  int dimension = 0;

  /// @brief Insert a value; return false if it is already in the span.
  bool insert(UInt value) {
    for (int bit = bits - 1; bit >= 0; bit--) {
      if (((value >> bit) & UInt(1)) == 0) {
        continue;
      }
      if (basis[bit] != 0) {
        value ^= basis[bit];
        continue;
      }
      basis[bit] = value;
      for (int lower = 0; lower < bit; lower++) {
        if ((basis[bit] >> lower) & UInt(1)) {
          basis[bit] ^= basis[lower];
        }
      }
      for (int higher = bit + 1; higher < bits; higher++) {
        if ((basis[higher] >> bit) & UInt(1)) {
          basis[higher] ^= basis[bit];
        }
      }
      dimension++;
      return true;
    }
    return false;
  }

  /// @brief Return whether value belongs to the represented xor span.
  bool contains(UInt value) const {
    for (int bit = bits - 1; bit >= 0; bit--) {
      if ((value >> bit) & UInt(1)) {
        value ^= basis[bit];
      }
    }
    return value == 0;
  }

  /// @brief Maximize seed xor x over all represented values x.
  UInt max_xor(UInt seed = 0) const {
    for (int bit = bits - 1; bit >= 0; bit--) {
      if ((seed ^ basis[bit]) > seed) {
        seed ^= basis[bit];
      }
    }
    return seed;
  }

  /// @brief Minimize seed xor x over all represented values x.
  UInt min_xor(UInt seed) const {
    for (int bit = bits - 1; bit >= 0; bit--) {
      if ((seed ^ basis[bit]) < seed) {
        seed ^= basis[bit];
      }
    }
    return seed;
  }

  /// @brief Return the k-th smallest distinct value in the span, including
  /// zero.
  UInt kth_smallest(std::uint64_t k) const {
    if (dimension < 64) {
      assert(k < (std::uint64_t(1) << dimension));
    }
    UInt result = 0;
    int index = 0;
    for (int bit = 0; bit < bits; bit++) {
      if (basis[bit] != 0) {
        if ((k >> index) & 1) {
          result ^= basis[bit];
        }
        index++;
      }
    }
    return result;
  }

  /// @brief Insert every basis vector from another span.
  void merge(const xor_basis &other) {
    for (UInt value : other.basis) {
      if (value != 0) {
        insert(value);
      }
    }
  }
};

/// @brief Return a basis of the intersection of two xor spans.  A homogeneous
/// system is built for U*a = V*b, with one equation per value bit. Gaussian
/// elimination produces a nullspace basis; projecting its U-coordinates gives
/// an independent basis of the intersection because both input lists are
/// required to be independent.
template <class UInt, int bits = std::numeric_limits<UInt>::digits>
std::vector<UInt>
xor_space_intersection(const std::vector<UInt> &first,
                       const std::vector<UInt> &second) {
  static_assert(std::is_integral_v<UInt> && std::is_unsigned_v<UInt>);
  static_assert(0 < bits && bits <= std::numeric_limits<UInt>::digits);
  const int left_size = int(first.size());
  const int variable_count = left_size + int(second.size());
  assert(variable_count <= 64);
  xor_basis<UInt, bits> left_check;
  xor_basis<UInt, bits> right_check;
  for (UInt value : first) {
    assert(left_check.insert(value));
  }
  for (UInt value : second) {
    assert(right_check.insert(value));
  }

  std::array<std::uint64_t, bits> equations{};
  for (int variable = 0; variable < variable_count; variable++) {
    UInt value = variable < left_size ? first[variable]
                                      : second[variable - left_size];
    for (int bit = 0; bit < bits; bit++) {
      if ((value >> bit) & UInt(1)) {
        equations[bit] |= std::uint64_t(1) << variable;
      }
    }
  }

  std::array<int, 64> pivot_row{};
  pivot_row.fill(-1);
  int rank = 0;
  for (int column = 0; column < variable_count; column++) {
    int pivot = rank;
    while (pivot < bits && ((equations[pivot] >> column) & 1) == 0) {
      pivot++;
    }
    if (pivot == bits) {
      continue;
    }
    std::swap(equations[rank], equations[pivot]);
    pivot_row[column] = rank;
    for (int row = 0; row < bits; row++) {
      if (row != rank && ((equations[row] >> column) & 1)) {
        equations[row] ^= equations[rank];
      }
    }
    rank++;
  }

  std::vector<UInt> result;
  for (int free_column = 0; free_column < variable_count; free_column++) {
    if (pivot_row[free_column] != -1) {
      continue;
    }
    std::uint64_t coefficients = std::uint64_t(1) << free_column;
    for (int column = 0; column < variable_count; column++) {
      int row = pivot_row[column];
      if (row != -1 && ((equations[row] >> free_column) & 1)) {
        coefficients |= std::uint64_t(1) << column;
      }
    }
    UInt value = 0;
    for (int index = 0; index < left_size; index++) {
      if ((coefficients >> index) & 1) {
        value ^= first[index];
      }
    }
    result.push_back(value);
  }
  return result;
}

} // namespace noya