Nyaan's Library

This documentation is automatically generated by online-judge-tools/verification-helper

View on GitHub

:heavy_check_mark: 任意mod二項係数
(modulo/arbitrary-mod-binomial.hpp)

任意mod二項係数

$\binom{n}{m} \mod M$を前計算$\mathrm{O}(\min(\frac{n\log M}{\log \log M},M))$、クエリ$\mathrm{O}(\frac{\log n \log M}{\log \log M})$で計算するライブラリ。

概要

論文　ferinさんの記事　hly1204さんの記事

$\binom{n}{m} \mod p^q$を計算出来ればCRTで$\binom{n}{m} \mod M$を復元できるので、ここでは素数冪に限定して議論する。

$\binom{n}{m} \pmod p$の公式

$\binom{n}{m} \mod p$は比較的容易に計算できる。

Lucasの定理

$\binom{n}{k} \mod p$は

\[n = a_r p^r + \cdots + a_1 p + a_0\] \[k = b_r p^r + \cdots + b_1 p + b_0\]

としたとき次の式で表される。

\[\binom{n}{k} \equiv \prod_{i=0}^r \binom{a_i}{b_i} \pmod p\]

(証明)

\[(1+x)^n \equiv \prod_{i=0}^r ((1+x)^{p^i})^{a_i}\] \[\equiv \prod_{i=0}^r (1 + x^{p^i})^{a_i} \equiv \prod_{i=0}^r \left(\sum_j \binom{a_i}{j} x^{jp^i}\right) \pmod p\]

と変形できるので、

\[\binom{n}{k} \equiv \lbrack x^k \rbrack (1 + x)^n\]

と合わせて上式を得る。

$\binom{n}{m}\mod p^k$の公式

法が素数冪の場合は様々な定理を用いて計算する必要がある。

Willsonの定理

素数$p$に対して次の合同式が成り立つ。

\[(p - 1)! \equiv -1 \pmod p\]

(証明) $\mod p$上で$1,\ldots,p-1$は逆元をただ一つ持つ。そのうち$a^2 \equiv 1 \pmod p \iff a \equiv \pm 1 \pmod p$を除いた$a=2,\ldots,p-2$は逆元を$2$から$p-2$の中に持つことから、$\prod_{a=2}^{p-2} a\equiv 1$を得る。よって$(p-1)! \equiv 1 \cdot (-1) \equiv -1 \pmod p$になる。

補題1

整数$n$に対して$(n!)_p$を$p$の倍数を除いた$n$以下の正整数の総積と置く。このとき素数冪$p^q$に対して

\[(p^q !)_p \equiv \delta \pmod {p^q}\]

が成り立つ。ただし$\delta$は$p=2,q\geq 3$の時は$1$、それ以外の時は$-1$とする。

(証明) Willsonの定理と同様に証明する。$\mod p^q$上で$p$と互いに素な数は唯一の逆元を持つため、$a^2 \not \equiv 1 \pmod {p^q}$である$a$の総積は$1$である。$a^2 \equiv 1 \pmod {p^q}$の解は

\[a=\left\lbrace \begin{array}{cc} 1 & p = 2 \wedge q = 1\newline \pm 1, 2^{q-1} \pm 1 & p = 2 \wedge q \geq 3 \newline \pm 1 & \mathrm{otherwise} \end{array} \right.\]

であるから、$\delta$は$p=2,q\geq 3$のとき$1$でそれ以外は$-1$になる。

また、$n \geq p^q$である$n$に対して補題を適用すると、$n = N_0 \bmod{p^q}$とおいて

\[(n!)_p \equiv \delta ^{\lfloor n / p^q \rfloor} (N_0!)_p \pmod{p ^ q}\]

という式が導かれる。

Legendreの定理ととKummerの定理

Legendreの定理とは次の有名な定理である。(証明略)

$n!$が$p$で割り切れる回数$\nu_p(n!)$は次の式で表される。

\[\nu_p(n!) =\sum_{1 \leq i} \left\lfloor \frac{n}{p^i} \right\rfloor\]

ここで次の補題2を考える。

$n$を$p$進数表示したときの各桁の和を$\sigma_p(n)$とおくと、$\sigma_p(n)$と$\nu_p(n!)$の間には次の関係式が成り立つ。

\[\nu_p(n!) = \frac{n - \sigma_p(n)}{p - 1}\]

証明は帰納的に行うことが出来る。

この補題を利用すると、次に説明するKummerの定理が証明できる。

$\binom{n}{m}$が$p$で割り切れる回数は$n-m$と$m$を$p$進数表示で足し算した時の繰り上がりの回数と等しい。

(証明)

$r = n - m$とおき、$p$進表示を$n = \sum_i n_i p^i$のように表すとする。また、$\epsilon_j$を$j$桁目で繰り上がりが起きた時$1$、起きないとき$0$を表す数とする。このとき、

\[n_j = m_j + r_j + \epsilon _{j-1} - p \epsilon_j\]

であり、

\[\nu_p\left(\frac{n!}{m!r!}\right) = \nu_p(n!) - \nu_p(m!) - \nu_p(r!)\] \[= \frac{\sigma_p(m) + \sigma_p(r) - \sigma_p(n)}{p - 1}\] \[= \sum_j \frac{m_j + r_j - n_j}{p - 1}\] \[= \sum_j \frac{\epsilon_j - p\epsilon_{j-1}}{p - 1} = \sum_j \epsilon_j\]

を得る。なお、$\epsilon_j$は次の式によって得られる。(Kummerの定理から従う。)

\[\epsilon_j = \left\lfloor \frac{n}{p^{j+1}}\right\rfloor - \left\lfloor \frac{m}{p^{j+1}}\right\rfloor - \left\lfloor \frac{r}{p^{j+1}}\right\rfloor\]

$n! \pmod {p^q}$の公式

正整数$n$と素数冪$p^q$に対して次の関係式が成り立つ。

\[{n!} / {p^{\sum_{j \geq 1} \lfloor n / p^j \rfloor }} \equiv \delta^{\sum_{j \geq q} \lfloor n / p^j \rfloor} \prod_{j \geq 0} (N_j !)_p \pmod{p^q}\]

ただし$N_j$は$\lfloor n / p^j \rfloor$と$\mod p^q$で合同な最小の非負整数とする。

(証明)

$j \geq 0$のとき

\[\lfloor n / p^j \rfloor ! / (p ^ {\lfloor n / p^{j + 1} \rfloor } \lfloor n / p^{j + 1} \rfloor ! ) \equiv (\lfloor n / p^{j} \rfloor !)_p\] \[\equiv \delta^{\lfloor n / p^{j + q} \rfloor} (N_j !)_p \pmod{p ^ q}\]

であるから、$j \geq 0$について上式を辺々掛ければ示せる。

$\binom{n}{m} \pmod {p^q}$の公式

上の式からただちに二項係数に対する以下の関係式が従う。

非負整数$n,m$と素数冪$p^q$について以下の式が成り立つ。

\[\frac{1}{p^{e_0}}\binom{n}{m} = \delta^{e_{q-1}} \prod_{j \geq 0} \frac{(N_j!)_p}{(M_j!)_p (R_j!)_p}\]

ただし$e_j$は$j$桁目以上で繰り上がりが発生した回数で、次の式で表される。

\[e_j = \sum_{j \lt i}\left(\left\lfloor \frac{n}{p^{i}}\right\rfloor - \left\lfloor \frac{m}{p^{i}}\right\rfloor - \left\lfloor \frac{r}{p^{i}}\right\rfloor\right)\]

$a \lt \min(n,p^q)$に対して$(a!)_p \bmod{p^q}$と$(a!)_p^{-1} \bmod{p^q}$を前計算しておくことで上の式は$\mathrm{O}(\log n)$で計算できる。

Barrett Reductionによる高速化

このライブラリではBarrett Reductionによる剰余算の高速化を行っているが、その正当性を示しておく。

参考: AtCoder Library

$n \lt B = 2^{64}, 1 \lt m \leq 2^{30}$のとき、

\[x = \left\lceil \frac{B}{m} \right\rceil\]

とおくと

\[\left\lfloor \frac{nx}{B} \right\rfloor = \left\lfloor \frac{n}{m} \right\rfloor + \varepsilon\]

が成り立つ。($\varepsilon$は$0$または$1$)

(証明) $x$の定義より

\[xm = B + r\ (0 \leq r \lt M)\]

が成り立つので、

\[0 \gt \frac{n}{m} - \frac{nx}{B} = n\left(\frac{1}{m} - \frac{x}{B} \right)\] \[=-\frac{n}{B}\cdot \frac{r}{M} \gt -1\]

であることから従う。(実用上は$m = 1$の時に$x$がオーバーフローすることに注意する必要がある。)

Depends on

Verified with

Code

#pragma once

#include <vector>

#include "../atcoder/math.hpp"
#include "../modint/barrett-reduction.hpp"

using namespace std;

#define PRIME_POWER_BINOMIAL_M_MAX ((1LL << 30) - 1)
#define PRIME_POWER_BINOMIAL_N_MAX 20000000

struct prime_power_binomial {
  int p, q, M;
  vector<int> fac, ifac, inv;
  int delta;
  Barrett bm, bp;

  prime_power_binomial(int _p, int _q) : p(_p), q(_q) {
    assert(1 < p && p <= PRIME_POWER_BINOMIAL_M_MAX);
    assert(_q > 0);
    long long m = 1;
    while (_q--) {
      m *= p;
      assert(m <= PRIME_POWER_BINOMIAL_M_MAX);
    }
    M = m;
    bm = Barrett(M), bp = Barrett(p);
    enumerate();
    delta = (p == 2 && q >= 3) ? 1 : M - 1;
  }

  void enumerate() {
    int MX = min<int>(M, PRIME_POWER_BINOMIAL_N_MAX + 10);
    fac.resize(MX);
    ifac.resize(MX);
    inv.resize(MX);
    fac[0] = ifac[0] = inv[0] = 1;
    fac[1] = ifac[1] = inv[1] = 1;
    for (int i = 2; i < MX; i++) {
      if (i % p == 0) {
        fac[i] = fac[i - 1];
        fac[i + 1] = bm.rem(1LL * fac[i - 1] * (i + 1));
        i++;
      } else {
        fac[i] = bm.rem(1LL * fac[i - 1] * i);
      }
    }
    ifac[MX - 1] = bm.pow(fac[MX - 1], M / p * (p - 1) - 1);
    for (int i = MX - 2; i > 1; --i) {
      if (i % p == 0) {
        ifac[i] = bm.rem(1LL * ifac[i + 1] * (i + 1));
        ifac[i - 1] = ifac[i];
        i--;
      } else {
        ifac[i] = bm.rem(1LL * ifac[i + 1] * (i + 1));
      }
    }
  }

  long long Lucas(long long n, long long m) {
    int res = 1;
    while (n) {
      int n0, m0;
      tie(n, n0) = bp.quorem(n);
      tie(m, m0) = bp.quorem(m);
      if (n0 < m0) return 0;
      res = bm.rem(1LL * res * fac[n0]);
      int buf = bm.rem(1LL * ifac[n0 - m0] * ifac[m0]);
      res = bm.rem(1LL * res * buf);
    }
    return res;
  }

  long long C(long long n, long long m) {
    if (n < m || n < 0 || m < 0) return 0;
    if (q == 1) return Lucas(n, m);
    long long r = n - m;
    int e0 = 0, eq = 0, i = 0;
    int res = 1;
    while (n) {
      res = bm.rem(1LL * res * fac[bm.rem(n)]);
      res = bm.rem(1LL * res * ifac[bm.rem(m)]);
      res = bm.rem(1LL * res * ifac[bm.rem(r)]);
      n = bp.quo(n);
      m = bp.quo(m);
      r = bp.quo(r);
      int eps = n - m - r;
      e0 += eps;
      if (e0 >= q) return 0;
      if (++i >= q) eq += eps;
    }
    if (eq & 1) res = bm.rem(1LL * res * delta);
    res = bm.rem(1LL * res * bm.pow(p, e0));
    return res;
  }
};

// constraints:
// (M <= 1e7 and max(N) <= 1e18) or (M < 2^30 and max(N) <= 2e7)
struct arbitrary_mod_binomial {
  int mod;
  vector<int> M;
  vector<prime_power_binomial> cs;

  arbitrary_mod_binomial(long long md) : mod(md) {
    assert(1 <= md);
    assert(md <= PRIME_POWER_BINOMIAL_M_MAX);
    for (int i = 2; i * i <= md; i++) {
      if (md % i == 0) {
        int j = 0, k = 1;
        while (md % i == 0) md /= i, j++, k *= i;
        M.push_back(k);
        cs.emplace_back(i, j);
        assert(M.back() == cs.back().M);
      }
    }
    if (md != 1) {
      M.push_back(md);
      cs.emplace_back(md, 1);
    }
    assert(M.size() == cs.size());
  }

  long long C(long long n, long long m) {
    if (mod == 1) return 0;
    vector<long long> rem, d;
    for (int i = 0; i < (int)cs.size(); i++) {
      rem.push_back(cs[i].C(n, m));
      d.push_back(M[i]);
    }
    return atcoder::crt(rem, d).first;
  }
};

#undef PRIME_POWER_BINOMIAL_M_MAX
#undef PRIME_POWER_BINOMIAL_N_MAX

/**
 * @brief 任意mod二項係数
 * @docs docs/modulo/arbitrary-mod-binomial.md
 */
#line 2 "modulo/arbitrary-mod-binomial.hpp"

#include <vector>

#line 1 "atcoder/math.hpp"



#include <algorithm>
#include <cassert>
#include <tuple>
#line 8 "atcoder/math.hpp"

#line 1 "atcoder/internal_math.hpp"



#include <utility>

#ifdef _MSC_VER
#include <intrin.h>
#endif

namespace atcoder {

namespace internal {

// @param m `1 <= m`
// @return x mod m
constexpr long long safe_mod(long long x, long long m) {
    x %= m;
    if (x < 0) x += m;
    return x;
}

// Fast modular multiplication by barrett reduction
// Reference: https://en.wikipedia.org/wiki/Barrett_reduction
// NOTE: reconsider after Ice Lake
struct barrett {
    unsigned int _m;
    unsigned long long im;

    // @param m `1 <= m < 2^31`
    barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1) {}

    // @return m
    unsigned int umod() const { return _m; }

    // @param a `0 <= a < m`
    // @param b `0 <= b < m`
    // @return `a * b % m`
    unsigned int mul(unsigned int a, unsigned int b) const {
        // [1] m = 1
        // a = b = im = 0, so okay

        // [2] m >= 2
        // im = ceil(2^64 / m)
        // -> im * m = 2^64 + r (0 <= r < m)
        // let z = a*b = c*m + d (0 <= c, d < m)
        // a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
        // c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2
        // ((ab * im) >> 64) == c or c + 1
        unsigned long long z = a;
        z *= b;
#ifdef _MSC_VER
        unsigned long long x;
        _umul128(z, im, &x);
#else
        unsigned long long x =
            (unsigned long long)(((unsigned __int128)(z)*im) >> 64);
#endif
        unsigned int v = (unsigned int)(z - x * _m);
        if (_m <= v) v += _m;
        return v;
    }
};

// @param n `0 <= n`
// @param m `1 <= m`
// @return `(x ** n) % m`
constexpr long long pow_mod_constexpr(long long x, long long n, int m) {
    if (m == 1) return 0;
    unsigned int _m = (unsigned int)(m);
    unsigned long long r = 1;
    unsigned long long y = safe_mod(x, m);
    while (n) {
        if (n & 1) r = (r * y) % _m;
        y = (y * y) % _m;
        n >>= 1;
    }
    return r;
}

// Reference:
// M. Forisek and J. Jancina,
// Fast Primality Testing for Integers That Fit into a Machine Word
// @param n `0 <= n`
constexpr bool is_prime_constexpr(int n) {
    if (n <= 1) return false;
    if (n == 2 || n == 7 || n == 61) return true;
    if (n % 2 == 0) return false;
    long long d = n - 1;
    while (d % 2 == 0) d /= 2;
    constexpr long long bases[3] = {2, 7, 61};
    for (long long a : bases) {
        long long t = d;
        long long y = pow_mod_constexpr(a, t, n);
        while (t != n - 1 && y != 1 && y != n - 1) {
            y = y * y % n;
            t <<= 1;
        }
        if (y != n - 1 && t % 2 == 0) {
            return false;
        }
    }
    return true;
}
template <int n> constexpr bool is_prime = is_prime_constexpr(n);

// @param b `1 <= b`
// @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g
constexpr std::pair<long long, long long> inv_gcd(long long a, long long b) {
    a = safe_mod(a, b);
    if (a == 0) return {b, 0};

    // Contracts:
    // [1] s - m0 * a = 0 (mod b)
    // [2] t - m1 * a = 0 (mod b)
    // [3] s * |m1| + t * |m0| <= b
    long long s = b, t = a;
    long long m0 = 0, m1 = 1;

    while (t) {
        long long u = s / t;
        s -= t * u;
        m0 -= m1 * u;  // |m1 * u| <= |m1| * s <= b

        // [3]:
        // (s - t * u) * |m1| + t * |m0 - m1 * u|
        // <= s * |m1| - t * u * |m1| + t * (|m0| + |m1| * u)
        // = s * |m1| + t * |m0| <= b

        auto tmp = s;
        s = t;
        t = tmp;
        tmp = m0;
        m0 = m1;
        m1 = tmp;
    }
    // by [3]: |m0| <= b/g
    // by g != b: |m0| < b/g
    if (m0 < 0) m0 += b / s;
    return {s, m0};
}

// Compile time primitive root
// @param m must be prime
// @return primitive root (and minimum in now)
constexpr int primitive_root_constexpr(int m) {
    if (m == 2) return 1;
    if (m == 167772161) return 3;
    if (m == 469762049) return 3;
    if (m == 754974721) return 11;
    if (m == 998244353) return 3;
    int divs[20] = {};
    divs[0] = 2;
    int cnt = 1;
    int x = (m - 1) / 2;
    while (x % 2 == 0) x /= 2;
    for (int i = 3; (long long)(i)*i <= x; i += 2) {
        if (x % i == 0) {
            divs[cnt++] = i;
            while (x % i == 0) {
                x /= i;
            }
        }
    }
    if (x > 1) {
        divs[cnt++] = x;
    }
    for (int g = 2;; g++) {
        bool ok = true;
        for (int i = 0; i < cnt; i++) {
            if (pow_mod_constexpr(g, (m - 1) / divs[i], m) == 1) {
                ok = false;
                break;
            }
        }
        if (ok) return g;
    }
}
template <int m> constexpr int primitive_root = primitive_root_constexpr(m);

}  // namespace internal

}  // namespace atcoder


#line 10 "atcoder/math.hpp"

namespace atcoder {

long long pow_mod(long long x, long long n, int m) {
  assert(0 <= n && 1 <= m);
  if (m == 1) return 0;
  internal::barrett bt((unsigned int)(m));
  unsigned int r = 1, y = (unsigned int)(internal::safe_mod(x, m));
  while (n) {
    if (n & 1) r = bt.mul(r, y);
    y = bt.mul(y, y);
    n >>= 1;
  }
  return r;
}

long long inv_mod(long long x, long long m) {
  assert(1 <= m);
  auto z = internal::inv_gcd(x, m);
  assert(z.first == 1);
  return z.second;
}

// (rem, mod)
std::pair<long long, long long> crt(const std::vector<long long>& r,
                                    const std::vector<long long>& m) {
  assert(r.size() == m.size());
  int n = int(r.size());
  // Contracts: 0 <= r0 < m0
  long long r0 = 0, m0 = 1;
  for (int i = 0; i < n; i++) {
    assert(1 <= m[i]);
    long long r1 = internal::safe_mod(r[i], m[i]), m1 = m[i];
    if (m0 < m1) {
      std::swap(r0, r1);
      std::swap(m0, m1);
    }
    if (m0 % m1 == 0) {
      if (r0 % m1 != r1) return {0, 0};
      continue;
    }
    // assume: m0 > m1, lcm(m0, m1) >= 2 * max(m0, m1)

    // (r0, m0), (r1, m1) -> (r2, m2 = lcm(m0, m1));
    // r2 % m0 = r0
    // r2 % m1 = r1
    // -> (r0 + x*m0) % m1 = r1
    // -> x*u0*g % (u1*g) = (r1 - r0) (u0*g = m0, u1*g = m1)
    // -> x = (r1 - r0) / g * inv(u0) (mod u1)

    // im = inv(u0) (mod u1) (0 <= im < u1)
    long long g, im;
    std::tie(g, im) = internal::inv_gcd(m0, m1);

    long long u1 = (m1 / g);
    // |r1 - r0| < (m0 + m1) <= lcm(m0, m1)
    if ((r1 - r0) % g) return {0, 0};

    // u1 * u1 <= m1 * m1 / g / g <= m0 * m1 / g = lcm(m0, m1)
    long long x = (r1 - r0) / g % u1 * im % u1;

    // |r0| + |m0 * x|
    // < m0 + m0 * (u1 - 1)
    // = m0 + m0 * m1 / g - m0
    // = lcm(m0, m1)
    r0 += x * m0;
    m0 *= u1;  // -> lcm(m0, m1)
    if (r0 < 0) r0 += m0;
  }
  return {r0, m0};
}

long long floor_sum(long long n, long long m, long long a, long long b) {
  long long ans = 0;
  if (a < 0) {
    unsigned long long a2 = internal::safe_mod(a, m);
    ans -= 1ULL * n * (n - 1) / 2 * ((a2 - a) / m);
    a = a2;
  }
  if (b < 0) {
    unsigned long long b2 = internal::safe_mod(b, m);
    ans -= 1ULL * n * ((b2 - b) / m);
    b = b2;
  }
  if (a >= m) {
    ans += (n - 1) * n * (a / m) / 2;
    a %= m;
  }
  if (b >= m) {
    ans += n * (b / m);
    b %= m;
  }
  long long y_max = (a * n + b) / m, x_max = (y_max * m - b);
  if (y_max == 0) return ans;
  ans += (n - (x_max + a - 1) / a) * y_max;
  ans += floor_sum(y_max, a, m, (a - x_max % a) % a);
  return ans;
}

}  // namespace atcoder


#line 2 "modint/barrett-reduction.hpp"

#line 4 "modint/barrett-reduction.hpp"
using namespace std;

struct Barrett {
  using u32 = unsigned int;
  using i64 = long long;
  using u64 = unsigned long long;
  u32 m;
  u64 im;
  Barrett() : m(), im() {}
  Barrett(int n) : m(n), im(u64(-1) / m + 1) {}
  constexpr inline i64 quo(u64 n) {
    u64 x = u64((__uint128_t(n) * im) >> 64);
    u32 r = n - x * m;
    return m <= r ? x - 1 : x;
  }
  constexpr inline i64 rem(u64 n) {
    u64 x = u64((__uint128_t(n) * im) >> 64);
    u32 r = n - x * m;
    return m <= r ? r + m : r;
  }
  constexpr inline pair<i64, int> quorem(u64 n) {
    u64 x = u64((__uint128_t(n) * im) >> 64);
    u32 r = n - x * m;
    if (m <= r) return {x - 1, r + m};
    return {x, r};
  }
  constexpr inline i64 pow(u64 n, i64 p) {
    u32 a = rem(n), r = m == 1 ? 0 : 1;
    while (p) {
      if (p & 1) r = rem(u64(r) * a);
      a = rem(u64(a) * a);
      p >>= 1;
    }
    return r;
  }
};
#line 7 "modulo/arbitrary-mod-binomial.hpp"

using namespace std;

#define PRIME_POWER_BINOMIAL_M_MAX ((1LL << 30) - 1)
#define PRIME_POWER_BINOMIAL_N_MAX 20000000

struct prime_power_binomial {
  int p, q, M;
  vector<int> fac, ifac, inv;
  int delta;
  Barrett bm, bp;

  prime_power_binomial(int _p, int _q) : p(_p), q(_q) {
    assert(1 < p && p <= PRIME_POWER_BINOMIAL_M_MAX);
    assert(_q > 0);
    long long m = 1;
    while (_q--) {
      m *= p;
      assert(m <= PRIME_POWER_BINOMIAL_M_MAX);
    }
    M = m;
    bm = Barrett(M), bp = Barrett(p);
    enumerate();
    delta = (p == 2 && q >= 3) ? 1 : M - 1;
  }

  void enumerate() {
    int MX = min<int>(M, PRIME_POWER_BINOMIAL_N_MAX + 10);
    fac.resize(MX);
    ifac.resize(MX);
    inv.resize(MX);
    fac[0] = ifac[0] = inv[0] = 1;
    fac[1] = ifac[1] = inv[1] = 1;
    for (int i = 2; i < MX; i++) {
      if (i % p == 0) {
        fac[i] = fac[i - 1];
        fac[i + 1] = bm.rem(1LL * fac[i - 1] * (i + 1));
        i++;
      } else {
        fac[i] = bm.rem(1LL * fac[i - 1] * i);
      }
    }
    ifac[MX - 1] = bm.pow(fac[MX - 1], M / p * (p - 1) - 1);
    for (int i = MX - 2; i > 1; --i) {
      if (i % p == 0) {
        ifac[i] = bm.rem(1LL * ifac[i + 1] * (i + 1));
        ifac[i - 1] = ifac[i];
        i--;
      } else {
        ifac[i] = bm.rem(1LL * ifac[i + 1] * (i + 1));
      }
    }
  }

  long long Lucas(long long n, long long m) {
    int res = 1;
    while (n) {
      int n0, m0;
      tie(n, n0) = bp.quorem(n);
      tie(m, m0) = bp.quorem(m);
      if (n0 < m0) return 0;
      res = bm.rem(1LL * res * fac[n0]);
      int buf = bm.rem(1LL * ifac[n0 - m0] * ifac[m0]);
      res = bm.rem(1LL * res * buf);
    }
    return res;
  }

  long long C(long long n, long long m) {
    if (n < m || n < 0 || m < 0) return 0;
    if (q == 1) return Lucas(n, m);
    long long r = n - m;
    int e0 = 0, eq = 0, i = 0;
    int res = 1;
    while (n) {
      res = bm.rem(1LL * res * fac[bm.rem(n)]);
      res = bm.rem(1LL * res * ifac[bm.rem(m)]);
      res = bm.rem(1LL * res * ifac[bm.rem(r)]);
      n = bp.quo(n);
      m = bp.quo(m);
      r = bp.quo(r);
      int eps = n - m - r;
      e0 += eps;
      if (e0 >= q) return 0;
      if (++i >= q) eq += eps;
    }
    if (eq & 1) res = bm.rem(1LL * res * delta);
    res = bm.rem(1LL * res * bm.pow(p, e0));
    return res;
  }
};

// constraints:
// (M <= 1e7 and max(N) <= 1e18) or (M < 2^30 and max(N) <= 2e7)
struct arbitrary_mod_binomial {
  int mod;
  vector<int> M;
  vector<prime_power_binomial> cs;

  arbitrary_mod_binomial(long long md) : mod(md) {
    assert(1 <= md);
    assert(md <= PRIME_POWER_BINOMIAL_M_MAX);
    for (int i = 2; i * i <= md; i++) {
      if (md % i == 0) {
        int j = 0, k = 1;
        while (md % i == 0) md /= i, j++, k *= i;
        M.push_back(k);
        cs.emplace_back(i, j);
        assert(M.back() == cs.back().M);
      }
    }
    if (md != 1) {
      M.push_back(md);
      cs.emplace_back(md, 1);
    }
    assert(M.size() == cs.size());
  }

  long long C(long long n, long long m) {
    if (mod == 1) return 0;
    vector<long long> rem, d;
    for (int i = 0; i < (int)cs.size(); i++) {
      rem.push_back(cs[i].C(n, m));
      d.push_back(M[i]);
    }
    return atcoder::crt(rem, d).first;
  }
};

#undef PRIME_POWER_BINOMIAL_M_MAX
#undef PRIME_POWER_BINOMIAL_N_MAX

/**
 * @brief 任意mod二項係数
 * @docs docs/modulo/arbitrary-mod-binomial.md
 */
Back to top page