Skip to content

montmort.hpp

SECTIONMath INCLUDEnoya/montmort.hpp

Return derangement counts D_0 through D_n modulo modulus. Separating the position occupied by the image of the first element gives D_n = (n-1)(D_{n-1}+D_{n-2}), with D_0=1 and D_1=0.

Verified by montmort_number_mod.

\[ \displaystyle D_n=(n-1)(D_{n-1}+D_{n-2}),\quad D_0=1,\ D_1=0 \]

Implementation

View on GitHub

#ifndef NOYA_MONTMORT_HPP
#define NOYA_MONTMORT_HPP 1

/// @complexity Time: O(n).  Space: O(n) for all values, or O(1) when only the
/// current recurrence state is retained.

#include <cassert>
#include <cstdint>
#include <vector>

namespace noya {

/// @brief Return derangement counts D_0 through D_n modulo modulus.
/// Separating the position occupied by the image of the first element gives
/// D_n = (n-1)(D_{n-1}+D_{n-2}), with D_0=1 and D_1=0.
inline std::vector<std::uint64_t> montmort_numbers(int n,
                                                  std::uint64_t modulus) {
  assert(n >= 0 && modulus >= 1);
  std::vector<std::uint64_t> result(n + 1);
  result[0] = 1 % modulus;
  for (int size = 2; size <= n; size++) {
    auto sum = (static_cast<unsigned __int128>(result[size - 1]) +
                result[size - 2]) % modulus;
    result[size] = std::uint64_t(
        static_cast<unsigned __int128>(size - 1) * sum % modulus);
  }
  return result;
}

} // namespace noya

#endif // NOYA_MONTMORT_HPP
#include <cassert>
#include <cstdint>
#include <vector>

/// @complexity Time: O(n).  Space: O(n) for all values, or O(1) when only the
/// current recurrence state is retained.

namespace noya {

/// @brief Return derangement counts D_0 through D_n modulo modulus.
/// Separating the position occupied by the image of the first element gives
/// D_n = (n-1)(D_{n-1}+D_{n-2}), with D_0=1 and D_1=0.
inline std::vector<std::uint64_t> montmort_numbers(int n,
                                                  std::uint64_t modulus) {
  assert(n >= 0 && modulus >= 1);
  std::vector<std::uint64_t> result(n + 1);
  result[0] = 1 % modulus;
  for (int size = 2; size <= n; size++) {
    auto sum = (static_cast<unsigned __int128>(result[size - 1]) +
                result[size - 2]) % modulus;
    result[size] = std::uint64_t(
        static_cast<unsigned __int128>(size - 1) * sum % modulus);
  }
  return result;
}

} // namespace noya