Quelle  bmat8.cpp   Sprache: C

//
// libsemigroups - C++ library for semigroups and monoids
// Copyright (C) 2019 Finn Smith
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
//

// This file contains an implemesntation of fast boolean matrices up to
// dimension 8 x 8.

#include "libsemigroups/bmat8.hpp"

#include <algorithm>  // for uniform_int_distribution
#include <array>      // for array
#include <climits>    // for CHAR_BIT
#include <iostream>   // for operator<<, ostringstream
#include <random>     // for mt19937, random_device
#include <vector>     // for vector

#include "libsemigroups/exception.hpp"  // for LIBSEMIGROUPS_EXCEPTION

namespace libsemigroups {

  ////////////////////////////////////////////////////////////////////////
  // Static helper data and functions
  ////////////////////////////////////////////////////////////////////////

  static std::array<uint64_t, 8> const ROW_MASK = {0xff00000000000000,
                                                   0xff000000000000,
                                                   0xff0000000000,
                                                   0xff00000000,
                                                   0xff000000,
                                                   0xff0000,
                                                   0xff00,
                                                   0xff};

  static std::array<uint64_t, 8> const COL_MASK = {0x8080808080808080,
                                                   0x4040404040404040,
                                                   0x2020202020202020,
                                                   0x1010101010101010,
                                                   0x808080808080808,
                                                   0x404040404040404,
                                                   0x202020202020202,
                                                   0x101010101010101};

  static std::array<uint64_t, 64> const BIT_MASK = {0x8000000000000000,
                                                    0x4000000000000000,
                                                    0x2000000000000000,
                                                    0x1000000000000000,
                                                    0x800000000000000,
                                                    0x400000000000000,
                                                    0x200000000000000,
                                                    0x100000000000000,
                                                    0x80000000000000,
                                                    0x40000000000000,
                                                    0x20000000000000,
                                                    0x10000000000000,
                                                    0x8000000000000,
                                                    0x4000000000000,
                                                    0x2000000000000,
                                                    0x1000000000000,
                                                    0x800000000000,
                                                    0x400000000000,
                                                    0x200000000000,
                                                    0x100000000000,
                                                    0x80000000000,
                                                    0x40000000000,
                                                    0x20000000000,
                                                    0x10000000000,
                                                    0x8000000000,
                                                    0x4000000000,
                                                    0x2000000000,
                                                    0x1000000000,
                                                    0x800000000,
                                                    0x400000000,
                                                    0x200000000,
                                                    0x100000000,
                                                    0x80000000,
                                                    0x40000000,
                                                    0x20000000,
                                                    0x10000000,
                                                    0x8000000,
                                                    0x4000000,
                                                    0x2000000,
                                                    0x1000000,
                                                    0x800000,
                                                    0x400000,
                                                    0x200000,
                                                    0x100000,
                                                    0x80000,
                                                    0x40000,
                                                    0x20000,
                                                    0x10000,
                                                    0x8000,
                                                    0x4000,
                                                    0x2000,
                                                    0x1000,
                                                    0x800,
                                                    0x400,
                                                    0x200,
                                                    0x100,
                                                    0x80,
                                                    0x40,
                                                    0x20,
                                                    0x10,
                                                    0x8,
                                                    0x4,
                                                    0x2,
                                                    0x1};

  // Cyclically shifts bits to left by 8m
  // https://stackoverflow.com/a/776523
  static uint64_t cyclic_shift(uint64_t n) noexcept {
    const unsigned int mask
        = (CHAR_BIT * sizeof(n) - 1);  // assumes width is a power of 2.
    // assert ( (c<=mask) &&"rotate by type width or more");
    unsigned int c = 8;
    c &= mask;
    return (n << c) | (n >> ((-c) & mask));
  }

  // Given 8 x 8 matrices A and B, this function returns a matrix where each
  // row is the corresponding row of B if that row of B is a subset of the
  // corresponding row of A, and 0 otherwise.
  static uint64_t zero_if_row_not_contained(uint64_t A, uint64_t B) noexcept {
    uint64_t w = A & B;
    for (size_t i = 0; i < 8; ++i) {
      if ((w & ROW_MASK[i]) == (B & ROW_MASK[i])) {
        w |= (B & ROW_MASK[i]);
      } else {
        w &= ~(ROW_MASK[i]);
      }
    }
    return w;
  }

  static std::array<uint64_t, 8> for_sorting = {0, 0, 0, 0, 0, 0, 0, 0};

  static inline void swap_for_sorting(size_t a, size_t b) noexcept {
    if (for_sorting[b] < for_sorting[a]) {
      std::swap(for_sorting[a], for_sorting[b]);
    }
  }

  ////////////////////////////////////////////////////////////////////////
  // BMat8 - static data members - private
  ////////////////////////////////////////////////////////////////////////

  std::random_device                      BMat8::_rd;
  std::mt19937                            BMat8::_gen(_rd());
  std::uniform_int_distribution<uint64_t> BMat8::_dist(0, 0xffffffffffffffff);

  ////////////////////////////////////////////////////////////////////////
  // BMat8 - constructors - public
  ////////////////////////////////////////////////////////////////////////

  BMat8::BMat8(std::vector<std::vector<bool>> const& mat) {
    if (0 == mat.size() || mat.size() > 8) {
      LIBSEMIGROUPS_EXCEPTION(
          "expected a vector with size in [1, 8], got a vector of size %d",
          mat.size());
    }
    _data        = 0;
    uint64_t pow = 1;
    pow          = pow << 63;
    for (auto row : mat) {
      if (row.size() != mat.size()) {
        LIBSEMIGROUPS_EXCEPTION("the vectors must all have the same length!");
      }
      for (auto entry : row) {
        if (entry) {
          _data ^= pow;
        }
        pow = pow >> 1;
      }
      pow = pow >> (8 - mat.size());
    }
  }

  // Not noexcept because it can throw.
  BMat8 BMat8::random(size_t dim) {
    if (0 == dim || dim > 8) {
      LIBSEMIGROUPS_EXCEPTION("the argument should be in [1, 8], got %d", dim);
    }
    BMat8 bm = BMat8::random();
    for (size_t i = dim; i < 8; ++i) {
      bm._data &= ~(ROW_MASK[i]);
      bm._data &= ~(COL_MASK[i]);
    }
    return bm;
  }

  ////////////////////////////////////////////////////////////////////////
  // BMat8 - member functions - public
  ////////////////////////////////////////////////////////////////////////

  BMat8 BMat8::row_space_basis() const noexcept {
    uint64_t out            = 0;
    uint64_t combined_masks = 0;

    BMat8 bm(_data);
    bm.sort_rows();
    uint64_t no_dups = bm._data;
    for (size_t i = 0; i < 7; ++i) {
      combined_masks |= ROW_MASK[i];
      while ((no_dups & ROW_MASK[i + 1]) << 8 == (no_dups & ROW_MASK[i])
             && (no_dups & ROW_MASK[i]) != 0) {
        no_dups = ((no_dups & combined_masks)
                   | (no_dups & ~combined_masks & ~ROW_MASK[i + 1]) << 8);
      }
    }
    uint64_t cm = no_dups;
    for (size_t i = 0; i < 7; ++i) {
      cm = cyclic_shift(cm);
      out |= zero_if_row_not_contained(no_dups, cm);
    }
    for (size_t i = 0; i < 8; ++i) {
      if ((out & ROW_MASK[i]) == (no_dups & ROW_MASK[i])) {
        out &= ~ROW_MASK[i];
      } else {
        out |= (no_dups & ROW_MASK[i]);
      }
    }
    combined_masks = 0;
    for (size_t i = 0; i < 8; ++i) {
      combined_masks |= ROW_MASK[i];
      while ((out & ROW_MASK[i]) == 0 && ((out & ~combined_masks) != 0)) {
        out = (out & combined_masks) | ((out & ~combined_masks) << 8);
      }
    }
    return BMat8(out);
  }

  void BMat8::set(size_t i, size_t j, bool val) {
    if (7 < i || 7 < j) {
      LIBSEMIGROUPS_EXCEPTION(
          "the arguments should each be at most 7, got %d as the %s argument",
          (7 < i) ? i : j,
          (7 < i) ? "first" : "second");
    }
    _data ^= (-val ^ _data) & BIT_MASK[8 * i + j];
  }

  BMat8 BMat8::operator*(BMat8 const& that) const noexcept {
    uint64_t y    = that.transpose()._data;
    uint64_t data = 0;
    uint64_t tmp  = 0;
    uint64_t diag = 0x8040201008040201;
    for (int i = 0; i < 8; ++i) {
      tmp = _data & y;
      tmp |= tmp >> 1;
      tmp |= tmp >> 2;
      tmp |= tmp >> 4;
      tmp &= 0x0101010101010101;
      tmp *= 255;
      tmp &= diag;
      data |= tmp;
      y    = cyclic_shift(y);
      diag = cyclic_shift(diag);
    }
    return BMat8(data);
  }

  size_t BMat8::row_space_size() const {
    std::array<char, 256> lookup;
    lookup.fill(0);
    std::vector<uint8_t> row_vec = row_space_basis().rows();
    auto                 last = std::remove(row_vec.begin(), row_vec.end(), 0);
    row_vec.erase(last, row_vec.end());
    for (uint8_t x : row_vec) {
      lookup[x] = true;
    }
    std::vector<uint8_t> row_space(row_vec.begin(), row_vec.end());
    for (size_t i = 0; i < row_space.size(); ++i) {
      for (uint8_t row : row_vec) {
        uint8_t x = row_space[i] | row;
        if (!lookup[x]) {
          row_space.push_back(x);
          lookup[x] = true;
        }
      }
    }
    return row_space.size() + 1;
  }

  std::vector<uint8_t> BMat8::rows() const {
    std::vector<uint8_t> rows;
    for (size_t i = 0; i < 8; ++i) {
      uint8_t row = static_cast<uint8_t>(_data << (8 * i) >> 56);
      rows.push_back(row);
    }
    return rows;
  }

  ////////////////////////////////////////////////////////////////////////
  // BMat8 - member functions - private
  ////////////////////////////////////////////////////////////////////////

  void BMat8::sort_rows() noexcept {
    for (size_t i = 0; i < 8; ++i) {
      for_sorting[i] = (_data << 8 * i) & ROW_MASK[0];
    }
    swap_for_sorting(0, 1);
    swap_for_sorting(2, 3);
    swap_for_sorting(0, 2);
    swap_for_sorting(1, 3);
    swap_for_sorting(1, 2);
    swap_for_sorting(4, 5);
    swap_for_sorting(6, 7);
    swap_for_sorting(4, 6);
    swap_for_sorting(5, 7);
    swap_for_sorting(5, 6);
    swap_for_sorting(0, 4);
    swap_for_sorting(1, 5);
    swap_for_sorting(1, 4);
    swap_for_sorting(2, 6);
    swap_for_sorting(3, 7);
    swap_for_sorting(3, 6);
    swap_for_sorting(2, 4);
    swap_for_sorting(3, 5);
    swap_for_sorting(3, 4);
    _data = 0;
    for (size_t i = 0; i < 7; ++i) {
      _data |= for_sorting[i];
      _data = (_data >> 8);
    }
    _data |= for_sorting[7];
  }

  ////////////////////////////////////////////////////////////////////////
  // BMat8 - friend functions
  ////////////////////////////////////////////////////////////////////////

  std::ostringstream& operator<<(std::ostringstream& os, BMat8 const& bm) {
    uint64_t x   = bm.to_int();
    uint64_t pow = 1;
    pow          = pow << 63;
    for (size_t i = 0; i < 8; ++i) {
      for (size_t j = 0; j < 8; ++j) {
        if (pow & x) {
          os << "1";
        } else {
          os << "0";
        }
        x = x << 1;
      }
      os << "\n";
    }
    return os;
  }

  ////////////////////////////////////////////////////////////////////////
  // BMat8 helpers
  ////////////////////////////////////////////////////////////////////////

  namespace bmat8_helpers {
    size_t minimum_dim(BMat8 const& x) noexcept {
      size_t   i = 0;
      uint64_t c = x.to_int();
      if (!c) {
        return 0;
      }
      uint64_t d = x.to_int();
      uint64_t y = x.transpose().to_int();
      uint64_t z = x.transpose().to_int();

      do {
        d = d >> 8;
        y = y >> 8;
        ++i;
      } while (i < 8 && (d << (8 * i)) == c && (y << (8 * i)) == z);

      return 9 - i;
    }
  }  // namespace bmat8_helpers
}  // namespace libsemigroups