LeviCivitaIterator.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <array>
8 #include <cstddef>
9 
10 #include "Utilities/Algorithm.hpp"
11 #include "Utilities/Array.hpp"
13 #include "Utilities/Numeric.hpp"
14 
15 /*!
16  * \ingroup DataStructuresGroup
17  *
18  * \brief Iterate over all nonzero index permutations for a Levi-Civita
19  * symbol.
20  *
21  * \details This class provides an iterator that allows you to loop over
22  * only the nonzero index permutations of a Levi-Civita symbol of dimension
23  * `dimension`. Inside the loop, the operator `()`
24  * returns an `std::array` containing an ordered list of the indices of this
25  * permutation, the operator `[]` returns a specific index from the same
26  * `std::array`, and the function `sign()` returns the sign of the
27  * Levi-Civita symbol for this permutation.
28  *
29  * \example
30  * \snippet Test_LeviCivitaIterator.cpp levi_civita_iterator_example
31  */
32 template <size_t Dim>
34  public:
35  explicit LeviCivitaIterator() noexcept = default;
36 
37  /// \cond HIDDEN_SYMBOLS
38  ~LeviCivitaIterator() = default;
39  // @{
40  /// No copy or move semantics
43  LeviCivitaIterator<Dim>& operator=(const LeviCivitaIterator<Dim>&) = delete;
44  LeviCivitaIterator<Dim>& operator=(LeviCivitaIterator<Dim>&&) = delete;
45  // @}
46  /// \endcond
47 
48  // Return false if the end of the loop is reached
49  explicit operator bool() const noexcept { return valid_; }
50 
51  // Increment the current permutation
52  LeviCivitaIterator& operator++() noexcept {
53  ++permutation_;
54  if (permutation_ >= signs_.size()) {
55  valid_ = false;
56  }
57  return *this;
58  }
59 
60  /// Return a `std::array` containing the current multi-index, an ordered list
61  /// of indices for the current permutation.
62  const std::array<size_t, Dim> operator()() const noexcept {
63  return static_cast<std::array<size_t, Dim>>(indexes_[permutation_]);
64  }
65 
66  /// Return a specific index from the multi-index of the current permutation.
67  const size_t& operator[](const size_t i) const noexcept {
68  return indexes_[permutation_][i];
69  }
70 
71  /// Return the sign of the Levi-Civita symbol for the current permutation.
72  int sign() const noexcept { return signs_[permutation_]; };
73 
74  private:
76  indexes() noexcept {
77  cpp17::array<cpp17::array<size_t, Dim>, factorial(Dim)> indexes{};
79  cpp2b::iota(index.begin(), index.end(), size_t(0));
80 
81  indexes[0] = index;
82 
83  // cpp20::next_permutation generates the different permuations of index
84  // loop over them to fill in the rest of the permutations in indexes
85  size_t permutation = 1;
86  while (cpp20::next_permutation(index.begin(), index.end())) {
87  indexes[permutation] = index;
88  ++permutation;
89  }
90  return indexes;
91  };
92 
93  static constexpr cpp17::array<int, factorial(Dim)> signs() noexcept {
95  // By construction, the sign of the first permutation is +1
96  signs[0] = 1;
97 
98  // How do you know whether the corresponding Levi Civita symbol is
99  // +1 or -1? To find out, compute the number, in the factoradic number
100  // system, corresponding to each permutation, and sum the digits. If the sum
101  // is even (odd), the corresponding Levi-Civita symbol will be +1 (-1). This
102  // works because there are Dim! unique permutations of the Levi Civita
103  // symbol indices, so each one can be represented by a
104  // (`Dim`)-digit factorial-number-system number. For more on the
105  // factoradic number system, see
106  // https://en.wikipedia.org/wiki/Factorial_number_system
107  auto factoradic_counter = cpp17::array<size_t, Dim>();
108  for (size_t permutation = 1; permutation < factorial(Dim); ++permutation) {
109  for (size_t i = 0; i < Dim; ++i) {
110  factoradic_counter[i]++;
111  if (factoradic_counter[i] < i + 1) {
112  break;
113  } else {
114  factoradic_counter[i] = 0;
115  }
116  }
117  signs[permutation] =
118  (cpp2b::accumulate(factoradic_counter.begin(),
119  factoradic_counter.end(), size_t(0)) %
120  2) == 0
121  ? 1
122  : -1;
123  }
124  return signs;
125  };
126 
127  // Note: here and throughout, use cpp17::array,
128  // which is constexpr (unlike std::array), to enable constexpr signs_
129  // and indexes_.
130  static constexpr cpp17::array<cpp17::array<size_t, Dim>, factorial(Dim)>
131  indexes_ = indexes();
132  static constexpr cpp17::array<int, factorial(Dim)> signs_{signs()};
133  size_t permutation_{0};
134  bool valid_{true};
135 };
constexpr void iota(ForwardIterator first, ForwardIterator last, T value)
Definition: Numeric.hpp:20
const size_t & operator[](const size_t i) const noexcept
Return a specific index from the multi-index of the current permutation.
Definition: LeviCivitaIterator.hpp:67
constexpr bool next_permutation(BidirectionalIterator first, BidirectionalIterator last, Compare comp)
Definition: Algorithm.hpp:93
const std::array< size_t, Dim > operator()() const noexcept
Return a std::array containing the current multi-index, an ordered list of indices for the current pe...
Definition: LeviCivitaIterator.hpp:62
Iterate over all nonzero index permutations for a Levi-Civita symbol.
Definition: LeviCivitaIterator.hpp:33
Define simple functions for constant expressions.
int sign() const noexcept
Return the sign of the Levi-Civita symbol for the current permutation.
Definition: LeviCivitaIterator.hpp:72
constexpr T accumulate(InputIt first, InputIt last, T init)
Definition: Numeric.hpp:33
constexpr uint64_t factorial(const uint64_t n) noexcept
Compute the factorial of .
Definition: ConstantExpressions.hpp:86
Definition: Array.hpp:26