|
1 | 1 | #ifndef SEQUANT_PERMUTATION_HPP |
2 | 2 | #define SEQUANT_PERMUTATION_HPP |
3 | 3 |
|
4 | | -#include <SeQuant/core/index.hpp> |
5 | 4 | #include <SeQuant/core/utility/macros.hpp> |
6 | 5 |
|
7 | 6 | #include <range/v3/algorithm.hpp> |
|
10 | 9 | #include <cstddef> |
11 | 10 | #include <cstdlib> |
12 | 11 | #include <set> |
13 | | -#include <type_traits> |
14 | | -#include <utility> |
| 12 | +#include <vector> |
15 | 13 |
|
16 | 14 | namespace sequant { |
17 | 15 |
|
18 | 16 | /// @brief Returns the number of cycles |
19 | 17 |
|
20 | 18 | /// Counts the number of cycles of a permutation represented in a 2-line form |
21 | 19 | /// by stacking \p v0 and \p v1 on top of each other. |
22 | | -/// @tparam Seq0 (reference to) a container type |
23 | | -/// @tparam Seq1 (reference to) a container type |
24 | | -/// @param v0 first sequence; if passed as an rvalue reference, it is moved from |
| 20 | +/// @tparam Seq0 a container type |
| 21 | +/// @tparam Seq1 a container type |
| 22 | +/// @param[in] v0 first sequence |
25 | 23 | /// @param[in] v1 second sequence |
26 | 24 | /// @pre \p v0 is a permutation of \p v1 |
27 | 25 | /// @return the number of cycles |
28 | 26 | template <typename Seq0, typename Seq1> |
29 | | -std::size_t count_cycles(Seq0&& v0, const Seq1& v1) { |
30 | | - std::remove_reference_t<Seq0> v(std::forward<Seq0>(v0)); |
31 | | - using T = std::decay_t<decltype(v[0])>; |
32 | | - SEQUANT_ASSERT(ranges::is_permutation(v, v1)); |
| 27 | +std::size_t count_cycles(const Seq0& v0, const Seq1& v1) { |
| 28 | + SEQUANT_ASSERT(ranges::is_permutation(v0, v1)); |
33 | 29 | // This function can't deal with duplicate entries in v0 or v1 |
34 | 30 | SEQUANT_ASSERT(std::set(std::begin(v0), std::end(v0)).size() == v0.size()); |
35 | 31 | SEQUANT_ASSERT(std::set(std::begin(v1), std::end(v1)).size() == v1.size()); |
36 | 32 |
|
37 | | - auto make_null = []() -> T { |
38 | | - if constexpr (std::is_arithmetic_v<T>) { |
39 | | - return -1; |
40 | | - } else if constexpr (std::is_same_v<T, Index>) { |
41 | | - return L"p_50"; |
42 | | - } |
43 | | - |
44 | | - SEQUANT_UNREACHABLE; |
45 | | - }; |
46 | | - |
47 | | - const auto null = make_null(); |
48 | | - SEQUANT_ASSERT(ranges::contains(v, null) == false); |
49 | | - SEQUANT_ASSERT(ranges::contains(v1, null) == false); |
50 | | - |
| 33 | + const auto n = v0.size(); |
| 34 | + std::vector<bool> visited(n, false); |
51 | 35 | std::size_t n_cycles = 0; |
52 | | - for (auto it = v.begin(); it != v.end(); ++it) { |
53 | | - if (*it != null) { |
54 | | - n_cycles++; |
55 | | - |
56 | | - auto idx = std::distance(v.begin(), it); |
57 | | - SEQUANT_ASSERT(idx >= 0); |
58 | | - |
59 | | - auto it0 = it; |
60 | | - |
61 | | - auto it1 = std::find(v1.begin(), v1.end(), *it0); |
62 | | - SEQUANT_ASSERT(it1 != v1.end()); |
63 | | - |
64 | | - auto idx1 = std::distance(v1.begin(), it1); |
65 | | - SEQUANT_ASSERT(idx1 >= 0); |
66 | 36 |
|
| 37 | + for (std::size_t i = 0; i < n; ++i) { |
| 38 | + if (!visited[i]) { |
| 39 | + ++n_cycles; |
| 40 | + auto j = i; |
67 | 41 | do { |
68 | | - it0 = std::find(v.begin(), v.end(), v[idx1]); |
69 | | - SEQUANT_ASSERT(it0 != v.end()); |
70 | | - |
71 | | - it1 = std::find(v1.begin(), v1.end(), *it0); |
72 | | - SEQUANT_ASSERT(it1 != v1.end()); |
73 | | - |
74 | | - idx1 = std::distance(v1.begin(), it1); |
75 | | - SEQUANT_ASSERT(idx1 >= 0); |
76 | | - |
77 | | - *it0 = null; |
78 | | - } while (idx1 != idx); |
| 42 | + visited[j] = true; |
| 43 | + auto it = std::find(std::begin(v1), std::end(v1), v0[j]); |
| 44 | + SEQUANT_ASSERT(it != std::end(v1)); |
| 45 | + j = static_cast<std::size_t>(std::distance(std::begin(v1), it)); |
| 46 | + } while (j != i); |
79 | 47 | } |
80 | 48 | } |
| 49 | + |
81 | 50 | return n_cycles; |
82 | | -}; |
| 51 | +} |
83 | 52 |
|
84 | 53 | /// computes parity of a permutation of 0 ... N-1 |
85 | 54 | /// |
|
0 commit comments