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

quality52%

¤ Dauer der Verarbeitung: 0.14 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.