順列の辞書順順位・復元

sequence/permutation_rank.hpp · GitHub で見る

共通ページ

0, ..., n-1 の順列について、辞書順で何番目かを求め、順位から元の順列を復元します。順位は 0 始まりです。入力は変更しません。

#include "lib/sequence/permutation_rank.hpp"
#include <cassert>

std::vector<int> p{2, 0, 3, 1};
auto rank = permutation_rank(p);
assert(rank == 13);
auto restored = permutation_unrank(4, rank);
assert(restored && *restored == p);
assert(!permutation_unrank(4, 24));

操作	引数・返り値
`permutation_rank<Rank = std::uint64_t>(p)`	順列 `p` の辞書順順位を `Rank` で返す
`permutation_unrank(n, rank)`	長さ `n` の順位 `rank` の順列を `std::optional<std::vector<int>>` で返す。整数型 `Rank` は引数から推論

p は std::span<const int> に変換できる vector<int>・array<int, n>・C 配列などを渡します。長さは INT_MAX 以下で、各整数が 0 以上 n 未満にちょうど 1 回ずつ現れるものとします。要素が範囲外・重複していれば std::invalid_argument、長さが上限を超えれば std::length_error を送出します。

順位は \(0\le r<n!\) です。空順列は順位 0 の 1 通りです。permutation_unrank は負の順位・n! 以上の順位には std::nullopt、n < 0 には std::invalid_argument を送出します。n! 自体が Rank に収まる必要はありません。

Rank には整数型を指定します。permutation_rank の結果が固定幅整数型に収まらない場合は std::overflow_error を送出します。既定型では長さ 20 まで全順列の順位が収まります。それ以上の長さでも、実際の順位が収まれば使えます。__int128_t・__uint128_t・boost::multiprecision::cpp_int や固定幅の多倍長整数も指定できます。順列の長さが Rank に収まる必要はありません。Boost は利用する側で include します。

#include "lib/sequence/permutation_rank.hpp"
#include <boost/multiprecision/cpp_int.hpp>
#include <cassert>
#include <numeric>

std::vector<int> p(50);
std::iota(p.rbegin(), p.rend(), 0);
auto rank = permutation_rank<boost::multiprecision::cpp_int>(p);
assert(permutation_unrank(50, rank) == p);

順位の剰余だけが必要なら、permutation_rank<Mint>(p) に modint 型を指定できます。型は mod() と整数からの構築・乗算代入・加算代入を備えるものとします。modint からの復元はできません。

Lehmer code と使用可能な値の個数を管理して計算します。階乗の前計算は不要です。

長さを \(n\) として、時間は \(O(n\log n)\)、追加空間は \(O(n)\) です。順位取得は \(n\le64\) なら \(O(n)\) 時間・\(O(1)\) 追加空間です。復元は Rank が 64 bit 以下の標準整数型なら \(O(n)\) 時間、返り値を除く追加空間は \(O(1)\) です。整数・modint の演算を \(O(1)\) とした計算量で、多倍長整数ではその演算コストが加わります。空順列は \(O(1)\) 時間・空間です。

コード

#line 1 "sequence/permutation_rank.hpp"

#include <array>
#include <bit>
#include <cstdint>
#include <limits>
#include <optional>
#include <span>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
#if defined(__BMI2__) && defined(__x86_64__)
#include <immintrin.h>
#endif

namespace permutation_order_detail {
inline constexpr auto factorial = [] {
	std::array<std::uint64_t, 21> result{};
	result[0] = 1;
	for (int i = 1; i <= 20; ++i) result[i] = result[i - 1] * i;
	return result;
}();
class Available {
	std::vector<int> tree;
  public:
	explicit Available(int n) : tree(std::size_t(n) + 1) {
		for (std::size_t i = 1; i < tree.size(); ++i) tree[i] = int(i & -i);
	}
	int prefix(int r) const {
		int result = 0;
		for (; r; r &= r - 1) result += tree[r];
		return result;
	}
	void erase(int p) {
		for (std::size_t i = std::size_t(p) + 1; i < tree.size(); i += i & -i) --tree[i];
	}
	int select(int k) const {
		std::size_t p = 0;
		for (std::size_t step = std::bit_floor(tree.size() - 1); step; step >>= 1)
			if (p + step < tree.size() && tree[p + step] <= k) {
				p += step;
				k -= tree[p];
			}
		return int(p);
	}
};

inline void restore(std::vector<int> &result) {
	int n = int(result.size()), prefix = 0;
	while (prefix < n && !result[prefix]) result[prefix] = prefix, ++prefix;
	std::uint64_t available = (std::uint64_t(1) << (n - prefix)) - 1;
	for (int &digit : std::span(result).subspan(prefix)) {
#if defined(__BMI2__) && defined(__x86_64__)
		digit = std::countr_zero(_pdep_u64(std::uint64_t(1) << digit, available));
#else
		std::uint64_t bits = available;
		for (int k = digit; k; --k) bits &= bits - 1;
		digit = std::countr_zero(bits);
#endif
		available &= ~(std::uint64_t(1) << digit);
		digit += prefix;
	}
}
}

template <class Rank = std::uint64_t>
Rank permutation_rank(std::span<const int> permutation) {
	static_assert(std::numeric_limits<Rank>::is_integer || requires { Rank::mod(); });
	if (permutation.size() > std::size_t(std::numeric_limits<int>::max())) throw std::length_error("permutation_rank");
	int n = int(permutation.size());
	Rank result = 0;
	auto append = [&](int digit, int base) {
		if constexpr (std::numeric_limits<Rank>::is_integer && std::numeric_limits<Rank>::is_bounded) {
			const Rank maximum = std::numeric_limits<Rank>::max();
			if constexpr (std::numeric_limits<Rank>::digits < std::numeric_limits<int>::digits) {
				if (digit > int(maximum) || (result != 0 && base > int(maximum))) throw std::overflow_error("permutation_rank");
				if (base > int(maximum)) { result = Rank(digit); return; }
			}
			if (result != 0 && result > (maximum - digit) / base) throw std::overflow_error("permutation_rank");
		}
		result *= base;
		result += digit;
	};
	if (n <= 64) {
		std::uint64_t remaining = n == 64 ? ~std::uint64_t(0) : (std::uint64_t(1) << n) - 1;
		std::uint64_t seen = 0;
		for (int x : permutation) {
			if (x < 0 || x >= n || (seen & (std::uint64_t(1) << x))) throw std::invalid_argument("permutation_rank");
			seen |= std::uint64_t(1) << x;
		}
		for (int i = 0; i < n; ++i) {
			int x = permutation[i];
			append(std::popcount(remaining & ((std::uint64_t(1) << x) - 1)), n - i);
			remaining ^= std::uint64_t(1) << x;
		}
	} else {
		permutation_order_detail::Available remaining(n);
		std::vector<bool> seen(n);
		for (int x : permutation) {
			if (x < 0 || x >= n || seen[x]) throw std::invalid_argument("permutation_rank");
			seen[x] = true;
		}
		for (int i = 0; i < n; ++i) {
			int x = permutation[i];
			append(remaining.prefix(x), n - i);
			remaining.erase(x);
		}
	}
	return result;
}

template <class Rank>
std::optional<std::vector<int>> permutation_unrank(int n, Rank rank) {
	static_assert(std::numeric_limits<Rank>::is_integer);
	if (n < 0) throw std::invalid_argument("permutation_unrank");
	if (rank < 0) return std::nullopt;
	if constexpr (std::is_integral_v<Rank> && sizeof(Rank) <= 8) {
		if (n <= 20 && std::uint64_t(rank) >= permutation_order_detail::factorial[n]) return std::nullopt;
		if (!n) return std::vector<int>{};
	}
	using Work = std::conditional_t<std::numeric_limits<Rank>::is_bounded &&
		(std::numeric_limits<Rank>::digits < std::numeric_limits<int>::digits), int, Rank>;
	Work remaining = static_cast<Work>(std::move(rank));
	std::vector<int> result(n);
	for (std::int64_t base = 2; base <= n && remaining != 0; ++base) {
		result[n - base] = int(remaining % base);
		remaining /= base;
	}
	if (remaining != 0) return std::nullopt;
	if constexpr (std::is_integral_v<Rank> && sizeof(Rank) <= 8) {
		permutation_order_detail::restore(result);
		return result;
	}
	permutation_order_detail::Available available(n);
	for (int &digit : result) {
		digit = available.select(digit);
		available.erase(digit);
	}
	return result;
}

#ifndef PERMUTATION_RANK_HPP
#define PERMUTATION_RANK_HPP

#include <array>
#include <bit>
#include <cstdint>
#include <limits>
#include <optional>
#include <span>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
#if defined(__BMI2__) && defined(__x86_64__)
#include <immintrin.h>
#endif

namespace permutation_order_detail {
inline constexpr auto factorial = [] {
	std::array<std::uint64_t, 21> result{};
	result[0] = 1;
	for (int i = 1; i <= 20; ++i) result[i] = result[i - 1] * i;
	return result;
}();
class Available {
	std::vector<int> tree;
  public:
	explicit Available(int n) : tree(std::size_t(n) + 1) {
		for (std::size_t i = 1; i < tree.size(); ++i) tree[i] = int(i & -i);
	}
	int prefix(int r) const {
		int result = 0;
		for (; r; r &= r - 1) result += tree[r];
		return result;
	}
	void erase(int p) {
		for (std::size_t i = std::size_t(p) + 1; i < tree.size(); i += i & -i) --tree[i];
	}
	int select(int k) const {
		std::size_t p = 0;
		for (std::size_t step = std::bit_floor(tree.size() - 1); step; step >>= 1)
			if (p + step < tree.size() && tree[p + step] <= k) {
				p += step;
				k -= tree[p];
			}
		return int(p);
	}
};

inline void restore(std::vector<int> &result) {
	int n = int(result.size()), prefix = 0;
	while (prefix < n && !result[prefix]) result[prefix] = prefix, ++prefix;
	std::uint64_t available = (std::uint64_t(1) << (n - prefix)) - 1;
	for (int &digit : std::span(result).subspan(prefix)) {
#if defined(__BMI2__) && defined(__x86_64__)
		digit = std::countr_zero(_pdep_u64(std::uint64_t(1) << digit, available));
#else
		std::uint64_t bits = available;
		for (int k = digit; k; --k) bits &= bits - 1;
		digit = std::countr_zero(bits);
#endif
		available &= ~(std::uint64_t(1) << digit);
		digit += prefix;
	}
}
}

template <class Rank = std::uint64_t>
Rank permutation_rank(std::span<const int> permutation) {
	static_assert(std::numeric_limits<Rank>::is_integer || requires { Rank::mod(); });
	if (permutation.size() > std::size_t(std::numeric_limits<int>::max())) throw std::length_error("permutation_rank");
	int n = int(permutation.size());
	Rank result = 0;
	auto append = [&](int digit, int base) {
		if constexpr (std::numeric_limits<Rank>::is_integer && std::numeric_limits<Rank>::is_bounded) {
			const Rank maximum = std::numeric_limits<Rank>::max();
			if constexpr (std::numeric_limits<Rank>::digits < std::numeric_limits<int>::digits) {
				if (digit > int(maximum) || (result != 0 && base > int(maximum))) throw std::overflow_error("permutation_rank");
				if (base > int(maximum)) { result = Rank(digit); return; }
			}
			if (result != 0 && result > (maximum - digit) / base) throw std::overflow_error("permutation_rank");
		}
		result *= base;
		result += digit;
	};
	if (n <= 64) {
		std::uint64_t remaining = n == 64 ? ~std::uint64_t(0) : (std::uint64_t(1) << n) - 1;
		std::uint64_t seen = 0;
		for (int x : permutation) {
			if (x < 0 || x >= n || (seen & (std::uint64_t(1) << x))) throw std::invalid_argument("permutation_rank");
			seen |= std::uint64_t(1) << x;
		}
		for (int i = 0; i < n; ++i) {
			int x = permutation[i];
			append(std::popcount(remaining & ((std::uint64_t(1) << x) - 1)), n - i);
			remaining ^= std::uint64_t(1) << x;
		}
	} else {
		permutation_order_detail::Available remaining(n);
		std::vector<bool> seen(n);
		for (int x : permutation) {
			if (x < 0 || x >= n || seen[x]) throw std::invalid_argument("permutation_rank");
			seen[x] = true;
		}
		for (int i = 0; i < n; ++i) {
			int x = permutation[i];
			append(remaining.prefix(x), n - i);
			remaining.erase(x);
		}
	}
	return result;
}

template <class Rank>
std::optional<std::vector<int>> permutation_unrank(int n, Rank rank) {
	static_assert(std::numeric_limits<Rank>::is_integer);
	if (n < 0) throw std::invalid_argument("permutation_unrank");
	if (rank < 0) return std::nullopt;
	if constexpr (std::is_integral_v<Rank> && sizeof(Rank) <= 8) {
		if (n <= 20 && std::uint64_t(rank) >= permutation_order_detail::factorial[n]) return std::nullopt;
		if (!n) return std::vector<int>{};
	}
	using Work = std::conditional_t<std::numeric_limits<Rank>::is_bounded &&
		(std::numeric_limits<Rank>::digits < std::numeric_limits<int>::digits), int, Rank>;
	Work remaining = static_cast<Work>(std::move(rank));
	std::vector<int> result(n);
	for (std::int64_t base = 2; base <= n && remaining != 0; ++base) {
		result[n - base] = int(remaining % base);
		remaining /= base;
	}
	if (remaining != 0) return std::nullopt;
	if constexpr (std::is_integral_v<Rank> && sizeof(Rank) <= 8) {
		permutation_order_detail::restore(result);
		return result;
	}
	permutation_order_detail::Available available(n);
	for (int &digit : result) {
		digit = available.select(digit);
		available.erase(digit);
	}
	return result;
}

#endif

#line 1 "sequence/permutation_rank.hpp"



#include <array>
#include <bit>
#include <cstdint>
#include <limits>
#include <optional>
#include <span>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
#if defined(__BMI2__) && defined(__x86_64__)
#include <immintrin.h>
#endif

namespace permutation_order_detail {
inline constexpr auto factorial = [] {
	std::array<std::uint64_t, 21> result{};
	result[0] = 1;
	for (int i = 1; i <= 20; ++i) result[i] = result[i - 1] * i;
	return result;
}();
class Available {
	std::vector<int> tree;
  public:
	explicit Available(int n) : tree(std::size_t(n) + 1) {
		for (std::size_t i = 1; i < tree.size(); ++i) tree[i] = int(i & -i);
	}
	int prefix(int r) const {
		int result = 0;
		for (; r; r &= r - 1) result += tree[r];
		return result;
	}
	void erase(int p) {
		for (std::size_t i = std::size_t(p) + 1; i < tree.size(); i += i & -i) --tree[i];
	}
	int select(int k) const {
		std::size_t p = 0;
		for (std::size_t step = std::bit_floor(tree.size() - 1); step; step >>= 1)
			if (p + step < tree.size() && tree[p + step] <= k) {
				p += step;
				k -= tree[p];
			}
		return int(p);
	}
};

inline void restore(std::vector<int> &result) {
	int n = int(result.size()), prefix = 0;
	while (prefix < n && !result[prefix]) result[prefix] = prefix, ++prefix;
	std::uint64_t available = (std::uint64_t(1) << (n - prefix)) - 1;
	for (int &digit : std::span(result).subspan(prefix)) {
#if defined(__BMI2__) && defined(__x86_64__)
		digit = std::countr_zero(_pdep_u64(std::uint64_t(1) << digit, available));
#else
		std::uint64_t bits = available;
		for (int k = digit; k; --k) bits &= bits - 1;
		digit = std::countr_zero(bits);
#endif
		available &= ~(std::uint64_t(1) << digit);
		digit += prefix;
	}
}
}

template <class Rank = std::uint64_t>
Rank permutation_rank(std::span<const int> permutation) {
	static_assert(std::numeric_limits<Rank>::is_integer || requires { Rank::mod(); });
	if (permutation.size() > std::size_t(std::numeric_limits<int>::max())) throw std::length_error("permutation_rank");
	int n = int(permutation.size());
	Rank result = 0;
	auto append = [&](int digit, int base) {
		if constexpr (std::numeric_limits<Rank>::is_integer && std::numeric_limits<Rank>::is_bounded) {
			const Rank maximum = std::numeric_limits<Rank>::max();
			if constexpr (std::numeric_limits<Rank>::digits < std::numeric_limits<int>::digits) {
				if (digit > int(maximum) || (result != 0 && base > int(maximum))) throw std::overflow_error("permutation_rank");
				if (base > int(maximum)) { result = Rank(digit); return; }
			}
			if (result != 0 && result > (maximum - digit) / base) throw std::overflow_error("permutation_rank");
		}
		result *= base;
		result += digit;
	};
	if (n <= 64) {
		std::uint64_t remaining = n == 64 ? ~std::uint64_t(0) : (std::uint64_t(1) << n) - 1;
		std::uint64_t seen = 0;
		for (int x : permutation) {
			if (x < 0 || x >= n || (seen & (std::uint64_t(1) << x))) throw std::invalid_argument("permutation_rank");
			seen |= std::uint64_t(1) << x;
		}
		for (int i = 0; i < n; ++i) {
			int x = permutation[i];
			append(std::popcount(remaining & ((std::uint64_t(1) << x) - 1)), n - i);
			remaining ^= std::uint64_t(1) << x;
		}
	} else {
		permutation_order_detail::Available remaining(n);
		std::vector<bool> seen(n);
		for (int x : permutation) {
			if (x < 0 || x >= n || seen[x]) throw std::invalid_argument("permutation_rank");
			seen[x] = true;
		}
		for (int i = 0; i < n; ++i) {
			int x = permutation[i];
			append(remaining.prefix(x), n - i);
			remaining.erase(x);
		}
	}
	return result;
}

template <class Rank>
std::optional<std::vector<int>> permutation_unrank(int n, Rank rank) {
	static_assert(std::numeric_limits<Rank>::is_integer);
	if (n < 0) throw std::invalid_argument("permutation_unrank");
	if (rank < 0) return std::nullopt;
	if constexpr (std::is_integral_v<Rank> && sizeof(Rank) <= 8) {
		if (n <= 20 && std::uint64_t(rank) >= permutation_order_detail::factorial[n]) return std::nullopt;
		if (!n) return std::vector<int>{};
	}
	using Work = std::conditional_t<std::numeric_limits<Rank>::is_bounded &&
		(std::numeric_limits<Rank>::digits < std::numeric_limits<int>::digits), int, Rank>;
	Work remaining = static_cast<Work>(std::move(rank));
	std::vector<int> result(n);
	for (std::int64_t base = 2; base <= n && remaining != 0; ++base) {
		result[n - base] = int(remaining % base);
		remaining /= base;
	}
	if (remaining != 0) return std::nullopt;
	if constexpr (std::is_integral_v<Rank> && sizeof(Rank) <= 8) {
		permutation_order_detail::restore(result);
		return result;
	}
	permutation_order_detail::Available available(n);
	for (int &digit : result) {
		digit = available.select(digit);
		available.erase(digit);
	}
	return result;
}