test/yosupo-judge/binomial_coefficient/main.test.cpp

• View this file on GitHub
• Last update: 2021-04-29 01:47:26+09:00
• Problem: https://judge.yosupo.jp/problem/binomial_coefficient
• Link: View error logs on GitHub Actions

Depends on

• Binomial coefficient (Mylib/Combinatorics/binomial_coefficient.cpp)
• Input tuple (Mylib/IO/input_tuple.cpp)
• Input tuples (Mylib/IO/input_tuples.cpp)
• Mod inverse (Mylib/Number/Mod/mod_inv.cpp)
• Mod pow (Mylib/Number/Mod/mod_pow.cpp)
• Chinese remainder theorem (Mylib/Number/chinese_remainder_algorithm.cpp)
• Extended Euclidean algorithm (Mylib/Number/extended_gcd.cpp)

Code

#define PROBLEM "https://judge.yosupo.jp/problem/binomial_coefficient"

#include "Mylib/Combinatorics/binomial_coefficient.cpp"
#include "Mylib/IO/input_tuples.cpp"

namespace hl = haar_lib;

int main() {
  std::cin.tie(0);
  std::ios::sync_with_stdio(false);

  int T, m;
  std::cin >> T >> m;

  hl::binomial_coefficient c(m);

  for (auto [n, k] : hl::input_tuples<int64_t, int64_t>(T)) {
    std::cout << c(n, k) << "\n";
  }

  return 0;
}

#line 1 "test/yosupo-judge/binomial_coefficient/main.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/binomial_coefficient"

#line 2 "Mylib/Number/Mod/mod_inv.cpp"
#include <cstdint>
#include <utility>

namespace haar_lib {
  constexpr int64_t mod_inv(int64_t a, int64_t m) {
    int64_t b = m, u = 1, v = 0;

    while (b) {
      int64_t t = a / b;
      a -= t * b;
      a = a ^ b;
      b = a ^ b;
      a = a ^ b;

      u -= t * v;
      u = u ^ v;
      v = u ^ v;
      u = u ^ v;
    }

    u %= m;
    if (u < 0) u += m;

    return u;
  }
}  // namespace haar_lib
#line 3 "Mylib/Number/Mod/mod_pow.cpp"

namespace haar_lib {
  constexpr int64_t mod_pow(int64_t n, int64_t p, int64_t m) {
    int64_t ret = 1;
    while (p > 0) {
      if (p & 1) (ret *= n) %= m;
      (n *= n) %= m;
      p >>= 1;
    }
    return ret;
  }
}  // namespace haar_lib
#line 2 "Mylib/Number/chinese_remainder_algorithm.cpp"
#include <cassert>
#include <optional>
#include <vector>
#line 2 "Mylib/Number/extended_gcd.cpp"
#include <tuple>

namespace haar_lib {
  auto ext_gcd(int64_t a, int64_t b) -> std::tuple<
      int64_t,  // gcd
      int64_t,  // p
      int64_t   // q
      > {
    if (b == 0) return std::make_tuple(a, 1, 0);
    const auto [d, q, p] = ext_gcd(b, (a + b) % b);
    return std::make_tuple(d, p, q - a / b * p);
  }
}  // namespace haar_lib
#line 6 "Mylib/Number/chinese_remainder_algorithm.cpp"

namespace haar_lib {
  std::optional<std::pair<int64_t, int64_t>> chinese_remainder_algorithm(
      int64_t b1, int64_t m1,
      int64_t b2, int64_t m2) {
    const auto [d, p, q] = ext_gcd(m1, m2);
    if ((b2 - b1) % d != 0) return std::nullopt;
    const int64_t m = m1 * m2 / d;
    const int64_t t = ((b2 - b1) * p / d) % (m2 / d);
    const int64_t r = (b1 + m1 * t + m) % m;
    return {{r, m}};
  }

  std::optional<std::pair<int64_t, int64_t>> chinese_remainder_algorithm(
      const std::vector<int64_t> &bs,
      const std::vector<int64_t> &ms) {
    assert(bs.size() == ms.size());
    int64_t R = 0, M = 1;
    for (int i = 0; i < (int) bs.size(); ++i) {
      const auto res = chinese_remainder_algorithm(R, M, bs[i], ms[i]);
      if (not res) return std::nullopt;
      const auto [r, m] = *res;
      R                 = r;
      M                 = m;
    }
    return {{R, M}};
  }
}  // namespace haar_lib
#line 6 "Mylib/Combinatorics/binomial_coefficient.cpp"

namespace haar_lib {
  class ext_lucas {
    std::vector<int64_t> prod_, inv_;
    int p_, q_;
    int64_t m_;

  public:
    ext_lucas(int p, int q) : p_(p), q_(q), m_(1) {
      for (int i = 0; i < q; ++i) m_ *= p_;

      prod_.assign(m_, 1);
      inv_.assign(m_, 1);

      for (int i = 1; i < m_; ++i) {
        prod_[i] = prod_[i - 1] * (i % p_ == 0 ? 1 : i) % m_;
      }

      inv_[m_ - 1] = mod_inv(prod_[m_ - 1], m_);
      for (int i = m_ - 1; i > 0; --i) {
        inv_[i - 1] = inv_[i] * (i % p_ == 0 ? 1 : i) % m_;
      }
    }

    int64_t operator()(int64_t n, int64_t k) const {
      int64_t r   = n - k;
      int64_t e   = 0;
      int64_t eq  = 0;
      int64_t ret = 1;

      for (int i = 0;;) {
        if (n == 0) { break; }

        (ret *= prod_[n % m_]) %= m_;
        (ret *= inv_[k % m_]) %= m_;
        (ret *= inv_[r % m_]) %= m_;

        n /= p_;
        k /= p_;
        r /= p_;

        e += n - k - r;

        if (e >= q_) return 0;

        i += 1;
        if (i >= q_) eq += n - k - r;
      }

      if ((p_ != 2 or q_ < 3) and eq % 2 == 1) ret = m_ - ret;

      (ret *= mod_pow(p_, e, m_)) %= m_;

      return ret;
    }
  };

  class binomial_coefficient {
    std::vector<std::pair<int, int>> m_primes;
    std::vector<ext_lucas> lu;
    std::vector<int64_t> ms;

  public:
    binomial_coefficient(int m) {
      for (int64_t i = 2LL; i * i <= m; ++i) {
        if (m % i == 0) {
          int t = 1, c = 0;
          while (m % i == 0) {
            m /= i;
            ++c;
            t *= i;
          }
          m_primes.emplace_back(i, c);
          ms.push_back(t);
        }
      }

      if (m != 1) {
        m_primes.emplace_back(m, 1);
        ms.push_back(m);
      }

      for (auto [p, q] : m_primes) {
        lu.push_back(ext_lucas(p, q));
      }
    }

    int64_t operator()(int64_t n, int64_t k) const {
      if (n < k or n < 0 or k < 0) return 0;

      std::vector<int64_t> bs;
      for (auto &a : lu) { bs.push_back(a(n, k)); }

      return chinese_remainder_algorithm(bs, ms).value().first;
    }
  };
}  // namespace haar_lib
#line 2 "Mylib/IO/input_tuples.cpp"
#include <initializer_list>
#include <iostream>
#line 6 "Mylib/IO/input_tuple.cpp"

namespace haar_lib {
  template <typename T, size_t... I>
  static void input_tuple_helper(std::istream &s, T &val, std::index_sequence<I...>) {
    (void) std::initializer_list<int>{(void(s >> std::get<I>(val)), 0)...};
  }

  template <typename T, typename U>
  std::istream &operator>>(std::istream &s, std::pair<T, U> &value) {
    s >> value.first >> value.second;
    return s;
  }

  template <typename... Args>
  std::istream &operator>>(std::istream &s, std::tuple<Args...> &value) {
    input_tuple_helper(s, value, std::make_index_sequence<sizeof...(Args)>());
    return s;
  }
}  // namespace haar_lib
#line 8 "Mylib/IO/input_tuples.cpp"

namespace haar_lib {
  template <typename... Args>
  class InputTuples {
    struct iter {
      using value_type = std::tuple<Args...>;
      value_type value;
      bool fetched = false;
      int N, c = 0;

      value_type operator*() {
        if (not fetched) {
          std::cin >> value;
        }
        return value;
      }

      void operator++() {
        ++c;
        fetched = false;
      }

      bool operator!=(iter &) const {
        return c < N;
      }

      iter(int N) : N(N) {}
    };

    int N;

  public:
    InputTuples(int N) : N(N) {}

    iter begin() const { return iter(N); }
    iter end() const { return iter(N); }
  };

  template <typename... Args>
  auto input_tuples(int N) {
    return InputTuples<Args...>(N);
  }
}  // namespace haar_lib
#line 5 "test/yosupo-judge/binomial_coefficient/main.test.cpp"

namespace hl = haar_lib;

int main() {
  std::cin.tie(0);
  std::ios::sync_with_stdio(false);

  int T, m;
  std::cin >> T >> m;

  hl::binomial_coefficient c(m);

  for (auto [n, k] : hl::input_tuples<int64_t, int64_t>(T)) {
    std::cout << c(n, k) << "\n";
  }

  return 0;
}

Back to top page