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             :  * \example
      30             :  * \snippet Test_LeviCivitaIterator.cpp levi_civita_iterator_example
      31             :  */
      32             : template <size_t Dim>
      33           1 : class LeviCivitaIterator {
      34             :  public:
      35           0 :   explicit LeviCivitaIterator() = default;
      36             : 
      37             :   /// \cond HIDDEN_SYMBOLS
      38             :   ~LeviCivitaIterator() = default;
      39             :   /// @{
      40             :   /// No copy or move semantics
      41             :   LeviCivitaIterator(const LeviCivitaIterator<Dim>&) = delete;
      42             :   LeviCivitaIterator(LeviCivitaIterator<Dim>&&) = delete;
      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           0 :   explicit operator bool() const { return valid_; }
      50             : 
      51             :   // Increment the current permutation
      52           0 :   LeviCivitaIterator& operator++() {
      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           1 :   const std::array<size_t, Dim> operator()() const {
      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           1 :   const size_t& operator[](const size_t i) const {
      68             :     return indexes_[permutation_][i];
      69             :   }
      70             : 
      71             :   /// Return the sign of the Levi-Civita symbol for the current permutation.
      72           1 :   int sign() const { return signs_[permutation_]; };
      73             : 
      74             :  private:
      75             :   static constexpr cpp20::array<cpp20::array<size_t, Dim>, factorial(Dim)>
      76           0 :   indexes() {
      77             :     cpp20::array<cpp20::array<size_t, Dim>, factorial(Dim)> indexes{};
      78             :     cpp20::array<size_t, Dim> index{};
      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           0 :   static constexpr cpp20::array<int, factorial(Dim)> signs() {
      94             :     cpp20::array<int, factorial(Dim)> signs{};
      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 = cpp20::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 cpp20::array, which is constexpr
     128             :   // (unlike std::array), to enable constexpr signs_ and indexes_.
     129             :   static constexpr cpp20::array<cpp20::array<size_t, Dim>, factorial(Dim)>
     130           0 :       indexes_ = indexes();
     131           0 :   static constexpr cpp20::array<int, factorial(Dim)> signs_{signs()};
     132           0 :   size_t permutation_{0};
     133           0 :   bool valid_{true};
     134             : };