Line data Source code
1 0 : // 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"
12 : #include "Utilities/ConstantExpressions.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 : * Note that the Levi-Civita symbol in odd dimensions (like 3) is invariant
30 : * under cyclic permutations of its indices, but in even dimensions (like 4)
31 : * it changes sign under cyclic permutations.
32 : *
33 : * \example
34 : * \snippet Test_LeviCivitaIterator.cpp levi_civita_iterator_example
35 : */
36 : template <size_t Dim>
37 1 : class LeviCivitaIterator {
38 : public:
39 0 : explicit LeviCivitaIterator() = default;
40 :
41 : /// \cond HIDDEN_SYMBOLS
42 : ~LeviCivitaIterator() = default;
43 : /// @{
44 : /// No copy or move semantics
45 : LeviCivitaIterator(const LeviCivitaIterator<Dim>&) = delete;
46 : LeviCivitaIterator(LeviCivitaIterator<Dim>&&) = delete;
47 : LeviCivitaIterator<Dim>& operator=(const LeviCivitaIterator<Dim>&) = delete;
48 : LeviCivitaIterator<Dim>& operator=(LeviCivitaIterator<Dim>&&) = delete;
49 : /// @}
50 : /// \endcond
51 :
52 : // Return false if the end of the loop is reached
53 0 : explicit operator bool() const { return valid_; }
54 :
55 : // Increment the current permutation
56 0 : LeviCivitaIterator& operator++() {
57 : ++permutation_;
58 : if (permutation_ >= signs_.size()) {
59 : valid_ = false;
60 : }
61 : return *this;
62 : }
63 :
64 : /// Return a `std::array` containing the current multi-index, an ordered list
65 : /// of indices for the current permutation.
66 1 : const std::array<size_t, Dim> operator()() const {
67 : return static_cast<std::array<size_t, Dim>>(indexes_[permutation_]);
68 : }
69 :
70 : /// Return a specific index from the multi-index of the current permutation.
71 1 : const size_t& operator[](const size_t i) const {
72 : return indexes_[permutation_][i];
73 : }
74 :
75 : /// Return the sign of the Levi-Civita symbol for the current permutation.
76 1 : int sign() const { return signs_[permutation_]; };
77 :
78 : private:
79 : static constexpr cpp20::array<cpp20::array<size_t, Dim>, factorial(Dim)>
80 0 : indexes() {
81 : cpp20::array<cpp20::array<size_t, Dim>, factorial(Dim)> indexes{};
82 : cpp20::array<size_t, Dim> index{};
83 : cpp2b::iota(index.begin(), index.end(), size_t(0));
84 :
85 : indexes[0] = index;
86 :
87 : // cpp20::next_permutation generates the different permutations of index
88 : // loop over them to fill in the rest of the permutations in indexes
89 : size_t permutation = 1;
90 : while (cpp20::next_permutation(index.begin(), index.end())) {
91 : indexes[permutation] = index;
92 : ++permutation;
93 : }
94 : return indexes;
95 : };
96 :
97 0 : static constexpr cpp20::array<int, factorial(Dim)> signs() {
98 : cpp20::array<int, factorial(Dim)> signs{};
99 : // By construction, the sign of the first permutation is +1
100 : signs[0] = 1;
101 :
102 : // How do you know whether the corresponding Levi Civita symbol is
103 : // +1 or -1? To find out, compute the number, in the factoradic number
104 : // system, corresponding to each permutation, and sum the digits. If the sum
105 : // is even (odd), the corresponding Levi-Civita symbol will be +1 (-1). This
106 : // works because there are Dim! unique permutations of the Levi Civita
107 : // symbol indices, so each one can be represented by a
108 : // (`Dim`)-digit factorial-number-system number. For more on the
109 : // factoradic number system, see
110 : // https://en.wikipedia.org/wiki/Factorial_number_system
111 : auto factoradic_counter = cpp20::array<size_t, Dim>();
112 : for (size_t permutation = 1; permutation < factorial(Dim); ++permutation) {
113 : for (size_t i = 0; i < Dim; ++i) {
114 : factoradic_counter[i]++;
115 : if (factoradic_counter[i] < i + 1) {
116 : break;
117 : } else {
118 : factoradic_counter[i] = 0;
119 : }
120 : }
121 : signs[permutation] =
122 : (cpp2b::accumulate(factoradic_counter.begin(),
123 : factoradic_counter.end(), size_t(0)) %
124 : 2) == 0
125 : ? 1
126 : : -1;
127 : }
128 : return signs;
129 : };
130 :
131 : // Note: here and throughout, use cpp20::array, which is constexpr
132 : // (unlike std::array), to enable constexpr signs_ and indexes_.
133 : static constexpr cpp20::array<cpp20::array<size_t, Dim>, factorial(Dim)>
134 0 : indexes_ = indexes();
135 0 : static constexpr cpp20::array<int, factorial(Dim)> signs_{signs()};
136 0 : size_t permutation_{0};
137 0 : bool valid_{true};
138 : };
|
