1 //****************************************************************************// 2 // Copyright (C) 2016 Florent Hivert <Florent.Hivert@lri.fr>, // 3 // // 4 // Distributed under the terms of the GNU General Public License (GPL) // 5 // // 6 // This code is distributed in the hope that it will be useful, // 7 // but WITHOUT ANY WARRANTY; without even the implied warranty of // 8 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU // 9 // General Public License for more details. // 10 // // 11 // The full text of the GPL is available at: // 12 // // 13 // http://www.gnu.org/licenses/ // 14 //****************************************************************************// 15 16 #ifndef HPCOMBI_PERM16_HPP_INCLUDED 17 #define HPCOMBI_PERM16_HPP_INCLUDED 18 19 #include <array> 20 #include <cassert> 21 #include <cstdint> 22 #include <functional> // less<> 23 #include <ostream> 24 #include <vector> 25 #include <x86intrin.h> 26 27 #include "epu.hpp" 28 #include "vect16.hpp" 29 30 namespace HPCombi { 31 32 // Forward declaration 33 struct Perm16; 34 struct PTransf16; 35 struct Transf16; 36 37 /** Partial transformation of @f$\{0\dots 15\}@f$ 38 * 39 */ 40 struct alignas(16) PTransf16 : public Vect16 { 41 sizeHPCombi::PTransf1642 static constexpr size_t size() { return 16; }; 43 44 using vect = HPCombi::Vect16; 45 using array = TPUBuild<epu8>::array; 46 47 PTransf16() = default; 48 HPCOMBI_CONSTEXPR_CONSTRUCTOR PTransf16(const PTransf16 &vv) = default; PTransf16HPCombi::PTransf1649 HPCOMBI_CONSTEXPR_CONSTRUCTOR PTransf16(const vect vv) : Vect16(vv) {} PTransf16HPCombi::PTransf1650 HPCOMBI_CONSTEXPR_CONSTRUCTOR PTransf16(const epu8 x) : Vect16(x) {} 51 PTransf16(std::vector<uint8_t> dom, std::vector<uint8_t> rng, 52 size_t = 0 /* unused */); 53 PTransf16(std::initializer_list<uint8_t> il); 54 55 PTransf16 &operator=(const PTransf16 &) = default; operator =HPCombi::PTransf1656 PTransf16 &operator=(const epu8 &vv) { v = vv; return *this; } 57 58 //! Return whether \c *this is a well constructed object validateHPCombi::PTransf1659 bool validate(size_t k = 16) const { 60 return HPCombi::is_partial_transformation(v, k); 61 } 62 63 //! The identity partial transformation. oneHPCombi::PTransf1664 static HPCOMBI_CONSTEXPR PTransf16 one() { return epu8id; } 65 //! The product of two partial transformations. operator *HPCombi::PTransf1666 PTransf16 operator*(const PTransf16 &p) const { 67 return HPCombi::permuted(v, p.v) | (p.v == Epu8(0xFF)); 68 } 69 70 /** Returns a mask for the image of \c *this */ 71 epu8 image_mask(bool complement=false) const; 72 /** Returns a bit mask for the image of \c *this */ 73 uint32_t image_bitset(bool complement=false) const; 74 /** Returns a mask for the domain of \c *this */ 75 epu8 domain_mask(bool complement=false) const; 76 /** Returns a bit mask for the domain of \c *this */ 77 uint32_t domain_bitset(bool complement=false) const; 78 79 /** Returns the partial right identity for \c *this */ 80 PTransf16 right_one() const; 81 /** Returns the partial left identity for \c *this */ 82 PTransf16 left_one() const; 83 84 /** Returns the size of the image of \c *this */ 85 uint32_t rank_ref() const; 86 /** Returns the size of the image of \c *this */ 87 uint32_t rank() const; 88 89 /** Returns a mask for the fix point of \c *this */ 90 epu8 fix_points_mask(bool complement=false) const; 91 /** Returns a bit mask for the fix point of \c *this */ 92 uint32_t fix_points_bitset(bool complement=false) const; 93 /** Returns the smallest fix point of \c *this */ 94 uint8_t smallest_fix_point() const; 95 /** Returns the smallest non fix point of \c *this */ 96 uint8_t smallest_moved_point() const; 97 /** Returns the largest fix point of \c *this */ 98 uint8_t largest_fix_point() const; 99 /** Returns the largest non fix point of \c *this */ 100 uint8_t largest_moved_point() const; 101 /** Returns the number of fix points of \c *this */ 102 uint8_t nb_fix_points() const; 103 }; 104 105 /** Full transformation of @f$\{0\dots 15\}@f$ 106 * 107 */ 108 struct Transf16 : public PTransf16 { 109 110 Transf16() = default; 111 HPCOMBI_CONSTEXPR_CONSTRUCTOR Transf16(const Transf16 &vv) = default; Transf16HPCombi::Transf16112 HPCOMBI_CONSTEXPR_CONSTRUCTOR Transf16(const vect vv) : PTransf16(vv) {} Transf16HPCombi::Transf16113 HPCOMBI_CONSTEXPR_CONSTRUCTOR Transf16(const epu8 x) : PTransf16(x) {} Transf16HPCombi::Transf16114 Transf16(std::initializer_list<uint8_t> il) : PTransf16(il) {} 115 Transf16 &operator=(const Transf16 &) = default; 116 117 //! Return whether \c *this is a well constructed object validateHPCombi::Transf16118 bool validate(size_t k = 16) const { 119 return HPCombi::is_transformation(v, k); 120 } 121 122 //! The identity transformation. oneHPCombi::Transf16123 static HPCOMBI_CONSTEXPR Transf16 one() { return epu8id; } 124 //! The product of two transformations. operator *HPCombi::Transf16125 Transf16 operator*(const Transf16 &p) const { 126 return HPCombi::permuted(v, p.v); 127 } 128 129 //! Construct a transformation from its 64 bits compressed. 130 explicit Transf16(uint64_t compressed); 131 //! The 64 bit compressed form of a transformation. 132 explicit operator uint64_t() const; 133 }; 134 135 /** Partial permutationof @f$\{0\dots 15\}@f$ 136 * 137 */ 138 struct PPerm16 : public PTransf16 { 139 140 PPerm16() = default; 141 HPCOMBI_CONSTEXPR_CONSTRUCTOR PPerm16(const PPerm16 &) = default; PPerm16HPCombi::PPerm16142 HPCOMBI_CONSTEXPR_CONSTRUCTOR PPerm16(const vect vv) : PTransf16(vv) {} PPerm16HPCombi::PPerm16143 HPCOMBI_CONSTEXPR_CONSTRUCTOR PPerm16(const epu8 x) : PTransf16(x) {} PPerm16HPCombi::PPerm16144 PPerm16(std::vector<uint8_t> dom, std::vector<uint8_t> rng, 145 size_t = 0 /* unused */) : PTransf16(dom, rng) {} PPerm16HPCombi::PPerm16146 PPerm16(std::initializer_list<uint8_t> il) : PTransf16(il) {} 147 PPerm16 &operator=(const PPerm16 &) = default; 148 149 //! Return whether \c *this is a well constructed object validateHPCombi::PPerm16150 bool validate(size_t k = 16) const { 151 return HPCombi::is_partial_permutation(v, k); 152 } 153 154 //! The identity partial permutations. oneHPCombi::PPerm16155 static HPCOMBI_CONSTEXPR PPerm16 one() { return epu8id; } 156 //! The product of two partial perrmutations. operator *HPCombi::PPerm16157 PPerm16 operator*(const PPerm16 &p) const { 158 return this->PTransf16::operator*(p); 159 } 160 161 /** @class common_inverse_pperm 162 * @brief The inverse of a partial permutation 163 * @details 164 * @returns the inverse of \c *this. The inverse of @f$p@f$ is the unique 165 * partial permutation @f$i@f$ such that @f$ p * i * p = p@f$ and 166 * @f$ i * p * i = i@f$ 167 * @par Example: 168 * @code 169 * Perm16 x = {0,3,2,4,0xFF,5,6,0xFF,8,9,11,0xFF,12,0xFF,0xFF,0xFF}; 170 * x.inverse() 171 * @endcode 172 * Returns 173 * @verbatim {0,0xFF,2,1,3,5,6,0xFF,8,9,0xFF,10,12,0xFF,0xFF,0xFF} @endverbatim 174 */ 175 /** @copydoc common_inverse_pperm 176 * @par Algorithm: 177 * @f$O(n)@f$ algorithm using reference cast to arrays 178 */ 179 PPerm16 inverse_ref() const; 180 /** @copydoc common_inverse_pperm 181 * @par Algorithm: 182 * @f$O(\log n)@f$ algorithm using some kind of vectorized dichotomic 183 * search. 184 */ 185 PPerm16 inverse_find() const; 186 right_oneHPCombi::PPerm16187 PPerm16 right_one() const { return PTransf16::right_one(); } left_oneHPCombi::PPerm16188 PPerm16 left_one() const { return PTransf16::left_one(); } 189 }; 190 191 /** Permutations of @f$\{0\dots 15\}@f$ 192 * 193 */ 194 struct Perm16 : public Transf16 /* public PPerm : diamond problem */ { 195 196 Perm16() = default; 197 HPCOMBI_CONSTEXPR_CONSTRUCTOR Perm16(const Perm16 &) = default; Perm16HPCombi::Perm16198 HPCOMBI_CONSTEXPR_CONSTRUCTOR Perm16(const vect vv) : Transf16(vv) {} Perm16HPCombi::Perm16199 HPCOMBI_CONSTEXPR_CONSTRUCTOR Perm16(const epu8 x) : Transf16(x) {} 200 Perm16 &operator=(const Perm16 &) = default; Perm16HPCombi::Perm16201 Perm16(std::initializer_list<uint8_t> il) : Transf16(il) {} 202 203 //! Return whether \c *this is a well constructed object validateHPCombi::Perm16204 bool validate(size_t k = 16) const { 205 return HPCombi::is_permutation(v, k); 206 } 207 208 // It's not possible to have a static constexpr member of same type as class 209 // being defined (see https://stackoverflow.com/questions/11928089/) 210 // therefore we chose to have functions. 211 //! The identity partial permutation. oneHPCombi::Perm16212 static HPCOMBI_CONSTEXPR Perm16 one() { return epu8id; } 213 //! The product of two permutations operator *HPCombi::Perm16214 Perm16 operator*(const Perm16 &p) const { 215 return HPCombi::permuted(v, p.v); 216 } 217 218 //! Construct a permutations from its 64 bits compressed. Perm16HPCombi::Perm16219 explicit Perm16(uint64_t compressed) : Transf16(compressed) {} 220 221 /** @class common_inverse 222 * @brief The inverse permutation 223 * @details 224 * @returns the inverse of \c *this 225 * @par Example: 226 * @code 227 * Perm16 x = {0,3,2,4,1,5,6,7,8,9,10,11,12,13,14,15}; 228 * x.inverse() 229 * @endcode 230 * Returns 231 * @verbatim {0,4,2,1,3,5,6,7,8,9,10,11,12,13,14,15} @endverbatim 232 */ 233 /** @copydoc common_inverse 234 * @par Algorithm: 235 * Reference @f$O(n)@f$ algorithm using loop and indexed access 236 */ 237 Perm16 inverse_ref() const; 238 /** @copydoc common_inverse 239 * @par Algorithm: 240 * @f$O(n)@f$ algorithm using reference cast to arrays 241 */ 242 Perm16 inverse_arr() const; 243 /** @copydoc common_inverse 244 * @par Algorithm: 245 * Insert the identity in the least significant bits and sort using a 246 * sorting network. The number of round of the optimal sorting network is 247 * open as far as I know, therefore the complexity is unknown. 248 */ 249 Perm16 inverse_sort() const; 250 /** @copydoc common_inverse 251 * @par Algorithm: 252 * @f$O(\log n)@f$ algorithm using some kind of vectorized dichotomic 253 * search. 254 */ inverse_findHPCombi::Perm16255 Perm16 inverse_find() const { return permutation_of(v, one()); } 256 /** @copydoc common_inverse 257 * @par Algorithm: 258 * 259 * Raise \e *this to power @f$\text{LCM}(1, 2, ..., n) - 1@f$ so complexity 260 * is in @f$O(log (\text{LCM}(1, 2, ..., n) - 1)) = O(n)@f$ 261 */ 262 Perm16 inverse_pow() const; 263 /** @copydoc common_inverse 264 * @par Algorithm: 265 * Compute power from @f$n/2@f$ to @f$n@f$, when @f$\sigma^k(i)=i@f$ then 266 * @f$\sigma^{-1}(i)=\sigma^{k-1}(i)@f$. Complexity @f$O(n)@f$ 267 */ 268 Perm16 inverse_cycl() const; 269 /** @copydoc common_inverse 270 * 271 * Frontend method: currently aliased to #inverse_cycl */ inverseHPCombi::Perm16272 Perm16 inverse() const { return inverse_cycl(); } 273 274 /** The elementary transposition exchanging @f$i@f$ and @f$i+1@f$ */ 275 static Perm16 elementary_transposition(uint64_t i); 276 /** A random permutation of size @f$n@f$*/ 277 static Perm16 random(uint64_t n = 16); 278 /** The \c r -th permutation of size \c n for the 279 * Steinhaus–Johnson–Trotter order. 280 */ 281 static Perm16 unrankSJT(int n, int r); 282 283 /** @class common_lehmer 284 * @brief The Lehmer code of a permutation 285 * @details 286 * @returns the Lehmer code of \c *this 287 * @par Example: 288 * @code 289 * Perm16 x = {0,3,2,4,1,5,6,7,8,9,10,11,12,13,14,15}; 290 * x.lehmer() 291 * @endcode 292 * Returns 293 * @verbatim {0,2,1,1,0,0,0,0,0,0,0,0,0,0,0,0} @endverbatim 294 */ 295 /** @copydoc common_lehmer 296 * @par Algorithm: 297 * Reference @f$O(n^2)@f$ algorithm using loop and indexed access 298 */ 299 epu8 lehmer_ref() const; 300 /** @copydoc common_lehmer 301 * @par Algorithm: 302 * Reference @f$O(n^2)@f$ algorithm using array, loop and indexed access 303 */ 304 epu8 lehmer_arr() const; 305 /** @copydoc common_lehmer 306 * @par Algorithm: 307 * Fast @f$O(n)@f$ algorithm using vector comparison 308 */ 309 epu8 lehmer() const; 310 311 /** @class common_length 312 * @brief The Coxeter length (ie: number of inversion) of a permutation 313 * @details 314 * @returns the number of inversions of \c *this 315 * @par Example: 316 * @code 317 * Perm16 x = {0,3,2,4,1,5,6,7,8,9,10,11,12,13,14,15}; 318 * x.length() 319 * @endcode 320 * Returns @verbatim 4 @endverbatim 321 */ 322 /** @copydoc common_length 323 * @par Algorithm: 324 * Reference @f$O(n^2)@f$ algorithm using loop and indexed access 325 */ 326 uint8_t length_ref() const; 327 /** @copydoc common_length 328 * @par Algorithm: 329 * Reference @f$O(n^2)@f$ algorithm using loop and indexed access after 330 * a cast to \c std::array 331 */ 332 uint8_t length_arr() const; 333 /** @copydoc common_length 334 * @par Algorithm: 335 * @f$O(n)@f$ using vector lehmer and fast horizontal sum 336 */ 337 uint8_t length() const; 338 339 /** @class common_nb_descent 340 * @brief The number of descent of a permutation 341 * @details 342 * @returns the number of inversions of \c *this 343 * @par Example: 344 * @code 345 * Perm16 x {0,3,2,4,1,5,6,7,8,9,10,11,12,13,14,15}; 346 * x.length() 347 * @endcode 348 * Returns @verbatim 2 @endverbatim 349 */ 350 /** @copydoc common_nb_descent 351 * @par Algorithm: 352 * Reference @f$O(n)@f$ using a loop 353 */ 354 uint8_t nb_descents_ref() const; 355 /** @copydoc common_nb_descent 356 * @par Algorithm: 357 * Reference @f$O(1)@f$ using vector shift and comparison 358 */ 359 uint8_t nb_descents() const; 360 361 /** The set partition of the cycles of a permutation 362 * @details 363 * @returns the a vector @f$v@f$ where @$fv[i]@$f contains the smallest 364 * element in the cycle of $i$ in \c *this 365 * @par Example: 366 * @code 367 * Perm16 x {1,2,3,6,0,5,4,7,8,9,10,11,12,15,14,13} 368 * x.cycles_partition() 369 * @endcode 370 * Returns 371 @verbatim 372 [ 0, 0, 0, 0, 0, 5, 0, 7, 8, 9,10,11,12,13,14,13] 373 @endverbatim 374 */ 375 epu8 cycles_partition() const; 376 377 /** @class common_nb_cycles 378 * @brief The number of cycles of a permutation 379 * @details 380 * @returns the number of cycles of \c *this 381 * @par Example: 382 * @code 383 * Perm16 x {1,2,3,6,0,5,4,7,8,9,10,11,12,15,14,13} 384 * x.nb_cycles() 385 * @endcode 386 * Returns @verbatim 10 @endverbatim 387 */ 388 /** @copydoc common_nb_cycles 389 * @par Algorithm: 390 * Reference @f$O(n)@f$ using a boolean vector 391 */ 392 uint8_t nb_cycles_ref() const; 393 /** @copydoc common_nb_cycles 394 * @par Algorithm: 395 * Reference @f$O(\log(n))@f$ using #cycles_partition 396 */ 397 uint8_t nb_cycles_unroll() const; 398 /** @copydoc common_nb_cycles 399 * @par Algorithm: aliased to #nb_cycles_unroll 400 */ nb_cyclesHPCombi::Perm16401 uint8_t nb_cycles() const { return nb_cycles_unroll(); } 402 403 /** @class common_left_weak_leq 404 * @brief Compare two permutations for the left weak order 405 * @par Example: 406 * @code 407 * Perm16 x{2,0,3,1}, y{3,0,2,1}; 408 * x.left_weak_leq(y) 409 * @endcode 410 * Returns @verbatim true @endverbatim 411 */ 412 /** @copydoc common_left_weak_leq 413 * @par Algorithm: 414 * Reference @f$O(n^2)@f$ testing inclusion of inversions one by one 415 */ 416 bool left_weak_leq_ref(Perm16 other) const; 417 /** @copydoc common_left_weak_leq 418 * @par Algorithm: 419 * Reference @f$O(n)@f$ with vectorized test of inclusion 420 */ 421 bool left_weak_leq_length(Perm16 other) const; 422 /** @copydoc common_left_weak_leq 423 * @par Algorithm: 424 * @f$O(n)@f$ algorithm using length 425 */ 426 bool left_weak_leq(Perm16 other) const; 427 428 }; 429 430 /*****************************************************************************/ 431 /** Memory layout concepts check ********************************************/ 432 /*****************************************************************************/ 433 434 static_assert(sizeof(epu8) == sizeof(Perm16), 435 "epu8 and Perm16 have a different memory layout !"); 436 static_assert(std::is_trivial<epu8>(), "epu8 is not a trivial class !"); 437 static_assert(std::is_trivial<Perm16>(), "Perm16 is not a trivial class !"); 438 439 } // namespace HPCombi 440 441 #include "perm16_impl.hpp" 442 443 namespace std { 444 445 template <> struct hash<HPCombi::PTransf16> { 446 //! A hash operator for #HPCombi::PTransf16 operator ()std::hash447 size_t operator()(const HPCombi::PTransf16 &ar) const { 448 return std::hash<HPCombi::epu8>{}(ar.v); 449 } 450 }; 451 452 template <> struct hash<HPCombi::Transf16> { 453 //! A hash operator for #HPCombi::Transf16 operator ()std::hash454 size_t operator()(const HPCombi::Transf16 &ar) const { 455 return uint64_t(ar); 456 } 457 }; 458 459 template <> struct hash<HPCombi::PPerm16> { 460 //! A hash operator for #HPCombi::PPerm16 operator ()std::hash461 size_t operator()(const HPCombi::PPerm16 &ar) const { 462 return std::hash<HPCombi::epu8>{}(ar.v); 463 } 464 }; 465 466 template <> struct hash<HPCombi::Perm16> { 467 //! A hash operator for #HPCombi::Perm16 operator ()std::hash468 size_t operator()(const HPCombi::Perm16 &ar) const { 469 return uint64_t(ar); 470 } 471 }; 472 473 } // namespace std 474 475 #endif // HPCOMBI_PERM16_HPP_INCLUDED 476