Ceizenpok’s Formula

题目描述

Dr. Ceizenp’ok from planet i1c5l became famous across the whole Universe thanks to his recent discovery – the Ceizenpok’s formula. This formula has only three arguments: $n, k$ and $m$, and its value is a number of $k$-combinations of a set of $n$ modulo $m$.
While the whole Universe is trying to guess what the formula is useful for, we need to automate its calculation.

题意概述

求${n \choose k} \bmod m$。
数据范围:$1 \le n \le 10^{18}, \; 0 \le k \le n, \; 2 \le m \le 10^6$。

算法分析


$$
x={n \choose k}, \; m=p_1^{a_1}p_2^{a_2} \cdots p_q^{a_q}
$$
可以列出方程组
$$
\left\{
\begin{array}{c}
x \equiv r_1 \pmod{p_1^{a_1}} \\
x \equiv r_2 \pmod{p_2^{a_2}} \\
\cdots \\
x \equiv r_q \pmod{p_q^{a_q}}
\end{array}
\right.
$$
由于模数两两互质,所以该方程组在模$m$意义下有唯一解。
考虑如何求$r_i$。实际上,我们要求的就是${n \choose k} \bmod p_i^{a_i}$。我们知道
$$
{n \choose k}={n! \over k!(n-k)!}
$$
那么只要求出$n! \bmod p_i^{a_i}, \; k! \bmod p_i^{a_i}, \; (n-k)! \bmod p_i^{a_i}$的值,就可以用逆元求出${n \choose k} \bmod p_i^{a_i}$。
对于如何求$n! \bmod p_i^{a_i}$,令
$$
f(n)=n! \bmod p_i^{a_i}
$$
由$x \equiv x+p_i^{a_i} \pmod{p_i^{a_i}}$,可得
$$
f(n)=\left(f\left(\left\lfloor{n \over p_i}\right\rfloor\right)\cdot p_i^{\lfloor n/p_i \rfloor}\cdot\left(\prod_{i \in [1, p_i^{a_i}], \; p_i\not\mid i}i\right)^{\lfloor n/p_i^{a_i} \rfloor}\cdot\prod_{i \in [1, n \bmod p_i^{a_i}], \; p_i\not\mid i}i\right)\bmod p_i^{a_i}
$$
但是$k!, \; (n-k)!$在模$p_i^{a_i}$意义下可能不存在逆元,因此需要将$n!, \; k!, \; (n-k)!$中的$p_i$因子提取出来,求出逆元后再乘回去。
这样就得到了所有$r_i$。用中国剩余定理求解方程组即可。

代码

/*
 * Your lucky number has been disconnected.
 */

#include <algorithm>
#include <cstdio>
#include <cstring>

#define int long long

int power(int a, int b, int m) {
  int ret = 1;
  for (; b; b >>= 1)
    b & 1 && ((ret *= a) %= m), (a *= a) %= m;
  return ret;
}

void ex_gcd(int a, int b, int &x, int &y) {
  if (!b) {
    x = 1, y = 0;
    return;
  }
  ex_gcd(b, a % b, y, x), y -= a / b * x;
}

int get_inv(int a, int b) {
  int x, y;
  ex_gcd(a, b, x, y);
  return (x + b) % b;
}

int get_fac(int n, int p, int k) {
  int m = power(p, k, 1e9), ret = 1;
  for (; n; n /= p) {
    if (n / m) {
      int rec = 1;
      for (int i = 2; i < m; ++i)
        if (i % p)
          (rec *= i) %= m;
      (ret *= power(rec, n / m, m)) %= m;
    }
    for (int i = n % m; i > 1; --i)
      if (i % p)
        (ret *= i) %= m;
  }
  return ret;
}

int calc(int N, int K, int M, int p, int k) {
  int a = get_fac(N, p, k), b = get_fac(K, p, k), c = get_fac(N - K, p, k),
      cnt = 0;
  for (int i = N; i; i /= p)
    cnt += i / p;
  for (int i = K; i; i /= p)
    cnt -= i / p;
  for (int i = N - K; i; i /= p)
    cnt -= i / p;
  int m = power(p, k, 1e9),
      ret = a * get_inv(b, m) % m * get_inv(c, m) % m * power(p, cnt, m) % m;
  return ret * (M / m) % M * get_inv(M / m, m) % M;
}

signed main() {
  int N, K, M, ans = 0;
  scanf("%lld%lld%lld", &N, &K, &M);
  for (int i = 2, t = M; t > 1; ++i)
    if (!(t % i)) {
      int k = 0;
      for (; !(t % i); ++k, t /= i)
        ;
      (ans += calc(N, K, M, i, k)) %= M;
    }
  printf("%lld\n", ans);
  return 0;
}

Notes on Module

逆元

对于模意义下的除法,$a/x \equiv b \pmod m \Leftrightarrow a \times x^{-1} \equiv b \pmod m$,则称$x^{-1}$为$x$在模$m$意义下的逆元。
关于逆元的求法可以参考Notes on Multiplicative Inverse
考虑线性同余方程$ax \equiv b \pmod m$,有
$$(a/(a, m))x \equiv b/(a, m) \pmod {m/(a, m)}$$
显然,当$(a, m) \not \mid b$时,原方程无解;否则,有
$$x \equiv (a/(a, m))^{-1} \times (b/(a, m)) \pmod {m/(a, m)}$$
因此原方程的解为
$$x \equiv (a/(a, m))^{-1} \times (b/(a, m)) + (m/(a, m)) \times k \pmod m \; (0 \le k \lt (a, m))$$
原方程可能有多解,也可能无解。


费马小定理

当$p$是素数时,对任意整数$x$都有$x^p \equiv x \pmod p$。若$p \not \mid x$,则
$$
\begin{align}
x^{p-1} &\equiv 1 \pmod p \\
x^{-1} &\equiv x^{p-2} \pmod p
\end{align}
$$
当$m$不是素数且$x$与$m$互质时,有
$$x^{\varphi(m)} \equiv 1 \pmod m$$
其中$\varphi(m)$为欧拉函数,其值等于不大于$m$的正整数中与$m$互质的数的个数。
利用线性筛,每次发现质数$p$时把它倍数的欧拉函数乘上$(p-1)/p$,就可以求出$1$到$n$的欧拉函数。


线性同余方程组

给定由若干个形如$a_ix \equiv b_i \pmod {m_i}$的方程组成的方程组,求它的解。如果有解,则一定有无穷多解,并可以表示为$x \equiv b \pmod m$的形式,因此问题转化为求解$b$和$m$。如果我们能求解方程组
$$
\left\{
\begin{align}
x &\equiv b_i \pmod {m_i} \\
ax &\equiv b_{i+1} \pmod {m_{i+1}}
\end{align}
\right.
$$
那么只要对方程逐个求解即可得到答案。
因为$x \equiv b_i \pmod {m_i}$,所以$x=tm_i+b_i$,代入第二个式子,得到
$$a(tm_i+b_i) \equiv b_{i+1} \pmod {m_{i+1}}$$
移项后得到
$$atm_i \equiv b_{i+1}-ab_i \pmod {m_{i+1}}$$
这是一个一次同余方程。当$(am_i, m_{i+1}) \not \mid b_{i+1}-ab_i$时原方程组无解。


中国剩余定理

如果同余方程组满足
$$\forall i \in [1, n], \; a_i=1 \land \forall i, j \in [1, n] \land i \neq j, \; (m_i, m_j)=1$$
那么答案的形式一定是
$$x \equiv b \pmod {\prod m_i}$$
反之,对于一个合数$n$,如果$n=ab \land (a, b)=1$,那么$x \bmod n$的值确定了之后,$x \bmod a$和$x \bmod b$的值也就确定了。因此,以$n$为模数来考虑和以$a$、$b$为模数来考虑是等价的。这个定理叫做中国剩余定理。
所以,对于模合数的情况只要考虑模若干$p^k$($p$为素数)的情况就可以了。如果$n$是一个square-free的数,那么问题就转化成模素数的情况,从而变得容易求解。
若$n$不是square-free的数,可以将每个数字$x$用一个二元组$(a, b)$表示,代表$x=ap_i^b$。有以下运算规则
$$
\begin{align}
(a, b)\times(c, d)&=(ac, b+d) \\
(a, b)/(c, d)&=(ac^{-1}, b-d)
\end{align}
$$
关于扩展中国剩余定理可参考Strange Way to Express Integers


$n!$

在计数问题中经常用到$n!$。因此有必要了解$n!$在模$p$意义下的一些性质。
假设$p$是素数,$n!=ap^e \; (a \not \mid p)$,求$a \bmod p$和$e$,$e$是$n!$能够迭代整除$p$的次数。显然$$e=n/p+n/p^2+n/p^3+ \cdots$$
因为$n/d$等于不超过$n$的能被$d$整除的正整数个数。
接下来计算$a \bmod p$。可以发现,不能被$p$整除的项在模$p$意义下呈周期性,它们的积为
$$((p-1)!)^{n/p} \times (n \bmod p)!$$
事实上,根据威尔逊定理,$(p-1)! \equiv -1 \pmod p$。因为除了$1$和$p-1$之外的其余项都可以和各自的逆元相乘得到$1$。
接下来计算可以被$p$整除的项的积。可以被$p$整除的项为$p, 2p, 3p, \ldots, (n/p)p$,将$p$消去后得到$1, 2, 3, \ldots, n/p$。因此问题的范围就由$n$缩小到了$n/p$。


${n \choose k} \bmod p$

当$p$不太大时,根据卢卡斯定理,令
$$n=\sum n_ip^i \; (0 \le n_i \lt p), \; k=\sum k_ip^i \; (0 \le k_i \lt p)$$

$${n \choose k} \equiv \prod {n_i \choose k_i} \pmod p$$
还有另一种思路。首先,把${n \choose k}$写成阶乘的形式
$${n \choose k}={n! \over k!(n-k)!}$$

$$n!=a_1p^{e_1}, \; k!=a_2p^{e_2}, \; (n-k)!=a_3p^{e_3}$$
从中可以看出,当$e_1 \gt e_2+e_3$时,${n \choose k}$可以被$p$整除;当$e_1=e_2+e_3$时,${n \choose k}$无法被$p$整除,此时${n \choose k} \equiv a_1(a_2a_3)^{-1} \pmod p$。
关于扩展卢卡斯定理可参考Ceizenpok’s Formula

Notes on Math Theory

Lucas

$n=(a_ka_{k-1}a_{k-2}\ldots a_0)_p, \; m=(b_kb_{k-1}b_{k-2}\ldots b_0)_p, \; {n \choose m} \equiv \prod_{i=0}^k {a_i \choose b_i} \pmod p$.
$n!$中$p$的次数$=\sum_{i=1}^{\infty} [{n \over p^i}]={n-S_p(n) \over p-1}$
${n \choose m}$中$p$的次数$p^\alpha \mid \mid {n \choose m} \Leftrightarrow p^\alpha \mid {n \choose m} \land p^{\alpha+1} \not \mid {n \choose m}$.
${n \choose m} \rightarrow p \not \mid {n \choose m} \Leftrightarrow {n-S_p(n) \over p-1}={m-S_p(m) \over p-1}+{n-m-S_(n-m) \over p-1} \Leftrightarrow m+(n-m)$在$p$进制下无进位。
$(1+x)^n \rightarrow x^m$系数${n \choose m}$.
$(1+x)^{p^k} \equiv 1+x^{p^k} \pmod p$.
$(1+x)^{a_kp^k+ \cdots +a_1p+a_0}=(1+x)^{a_kp^k} \cdots (1+x)^{a_1p}(1+x)^{a_0}$.
$\equiv (1+x^{p^k})^{a_k} \cdots (1+x^p)^{a_1}(1+x)^{a_0} \pmod p$.
$n!!=1 \times 3 \times 5 \times \cdots \times n$.


证明${n \choose k} (0 \le k \le n)$为奇数的个数是$2$的幂次。


Euler

$x^{\varphi(n)} \equiv 1 \pmod n, (x, n)=1$.
$x^n \equiv x^{n-\varphi(n)} \pmod n$.
$\bmod n$缩系.
$a_1 \times a_2 \times \cdots \times a_{\varphi(n)} \equiv a_1x \times a_2x \times \cdots \times a_{\varphi(n)}x \pmod n$.
$n$为素数,$x^{p-1} \equiv 1 \pmod p, p \not \mid x$.
$x^p \equiv x \pmod p$.
$\sum_{d \mid n} \varphi(d)=n$.
$\varphi(p^{\alpha})=p^{\alpha}-p^{\alpha-1}$.
$n=p_1^{\alpha_1} \cdots p_k^{\alpha_k}$.
$\sum_{\beta_i \le \alpha_i} \varphi(p_1^{\beta_1} \cdots p_k^{\beta_k})=\sum_{\beta_i \le \alpha_i} \prod_{i=1}^k \varphi(p_i^{\beta_i})=\prod_{i=1}^k \sum_{\beta_i \le \alpha_i} \varphi(p_i^{\beta_i})=\prod p_i^{\alpha_i}$.


求满足$n \mid 2^n-1$的$n$的个数。
$n=1$成立。
$n>1$时,设$p$为$n$最小的素因子。
满足$p \mid 2^{d_0}-1$的最小的$d_0$,$\forall p \mid 2^n-1, \; d_0 \mid n$.
$\begin{cases} p \mid 2^n-1\\ p \mid 2^{p-1}-1\\ p \mid 2^d-1 \end{cases} \Rightarrow \begin{cases} d \mid n\\ d \mid p-1 \end{cases} \Rightarrow d=1 \Rightarrow p \mid 1 \Rightarrow n=1$.
$ax+by=(a, b)$.
$nx+(p-1)y=(n, p-1)$.
$p \mid 2^{(n, p-1)}-1$.
∵$(n, p-1)=1$.
∴$p \mid 1$.


$1, x, x^2, \ldots, x^{\varphi(n)-1}$构成$\bmod n$的缩系,则称$x$为$n$的原根,记作$g$。
$p$是素数时,${g, g^2, \ldots, g^{p-1}} \equiv {1, 2, \ldots, p-1} \pmod p$.
$n$的原根有$\varphi(\varphi(n))$个。


给定$k$,求$\sum_{i=1}^{p-1} i^k \pmod p$。
${1, 2, \ldots, p-1} \Rightarrow {a, 2a, \ldots, (p-1)a}$.
$\sum i^k \equiv \sum (a_i)^k \Rightarrow (a^k-1) \sum i^k \equiv 0 \pmod p$.


二项同余式

$x^k \equiv n \pmod p$无解或有$(k, p-1)$个解。
$f(x) \equiv 0 \pmod p$.
$x^k \equiv 1 \pmod 1$有$d$个解,$d=(k, p-1)=ak+b(p-1)$。
$\begin{cases} x^{p-1} \equiv 1 \pmod p\\ x^k \equiv 1 \pmod 1 \end{cases} \Rightarrow x^{(k, p-1)} \equiv 1 \pmod p$至多$d$个解。
${x^{p-1}-1 \over x^d-1}=x^{p-1-d}+x^{p-1-2d}+ \cdots +x^d+1 \equiv 1 \pmod p$至多$p-1-d$个解。
$x^k \equiv 1 \pmod 1$恰有$d$个解。
$x \rightarrow p$的缩系,$x^k$有${p-1 \over (k, p-1)}$个不同的取值。
$x^{p-1} \equiv 1 \pmod p$且$x^d \not \equiv 1 \pmod p, d \lt p-1$,$x$是原根。
令$\gamma(n)$表示次数为$n$的数的个数。$p-1=q_1^{\alpha_1} \cdots q_k^{\alpha_k}, \gamma(q^{\alpha})=q^{\alpha}-q^{\alpha-1}$,且$\gamma$是积性函数。
$x^{q^{\alpha}} \equiv 1 \pmod p$.
$x^{q^{\alpha-1}} \not \equiv 1 \Leftrightarrow \forall d \lt q^{alpha}, \; x^d \not \equiv 1 \pmod p$.
$x^{(q^{\alpha}, d)} \equiv 1 \Leftrightarrow x^{q^{\beta}} \equiv 1 \pmod p$.
$\gamma(l_1l_2)=\gamma(l_1)\gamma(l_2)$.
$h_1$次数是$l_1$,$h_2$次数是$l_2 \Rightarrow h_1h_2$的次数是$l_1l_2$。
$x$次数是$l_1l_2 \Rightarrow x_1$次数是$l_1$,$x_2$次数是$l_2$。


Prime

$\lim \sum_{p=1}^{\infty} (1-{1 \over p})=0$.


$f(x) \equiv 0 \pmod {p^{\alpha}}$的解数$=f(x) \equiv 0 \pmod p$的解数$\Leftarrow \begin{cases} f(x) \equiv 0\\ f'(x) \equiv 0 \end{cases} \pmod p$无解。
$x^2+1 \equiv 0 \pmod {p^{\alpha}}, p \ge 3$无解或恰有$2$解。


勾股数的通式

$a^2+b^2=c^2 (a, b, c \gt 0)$.
$(a, b, c)=1$本原。
$a^2+b^2=c^2 \Rightarrow b^2=c^2-a^2=(c-a)(c+a)$.
$(c-a, c+a)=(2c, c-a)=(2, c-a)=1$.
$2 \not \mid c-a$.
设$b$为偶数,$a$为奇数,则$c$为奇数。
$\begin{cases} {c-a \over 2}=m^2\\ {c+a \over 2}=n^2 \end{cases}$.
$\begin{cases} b=2mn\\ a=n^2-m^2\\ c=n^2+m^2 \end{cases}$.


多项式

拉格朗日差值公式:$f(x)=\sum_{i=1}^n \prod_{j \neq i} {(x-x_j) \over (x_i-x_j)}f(x_i)$.
牛顿差值公式:$f(x)=a_0+a_1(x-x_1)+a_2(x-x_1)(x-x_2)+ \cdots +a_n \prod_{i=1}^n (x-x_i)$.


有理根 Eisenstein

$f \in Z[x], f(x)=a_nx^n+ \cdots +a_1x+a_0$.
$\pm {p \over q} \Rightarrow p \mid a_n, q \mid a_0$.
$f(x)=c_nx^n+c_{n-1}x^{n-1}+ \cdots +c_1x+c_0$可约$\Leftarrow p \mid c_0, p \mid c_1, \ldots, p \mid c_{n-1}, p \not \mid c_n, p^2 \not \mid c_0$.
$f(x)=c_nx^n+c_{n-1}x^{n-1}+ \cdots +c_1x+c_0$可约$\Leftarrow p \mid (c_0, c_1, \ldots, c_{n-1}), p \not \mid c_n, p^2 \not \mid c_0$.
$f(x)=x^{p-1}+x^{p-2}+ \cdots +x+1={x^p-1 \over {x-1}}$.
令$y=x-1$,$f(x)=g(y)={(y+1)^p-1 \over y}=y^{p-1}+{p \choose 1}y^{p-2}+ \cdots +{p \choose p-1}$.


升幂定理

令$V_p(n)$表示$n$中$p$的次数。
$5^{210000}-2^{210000}$里有多少$3$的幂次。$2$次。
$(1+p)^k-1$里有多少$p$的幂次。$1+V_p(k)$次。
$5^{33333}+1$里有多少$3$的幂次。$2$次。
$p^{\alpha} \mid \mid a-b, (a, p)=1, (b, p)=1, p \ge 3$.
$V_p(a^m-b^m)=\alpha+V_p(m)$.
设$p^{\alpha} \mid \mid a^n-b^n$.
$p^{\alpha} \mid \mid a^{pn}-b^{pn}$.
$p^{\alpha} \mid \mid a^{kn}-b^{kn}, (k, p)=1$.


二次剩余

$ax^2+bx+c \equiv 0 \pmod p$.
$x^2 \equiv a \pmod p$.
勒让德符号:$({a \over p})=\begin{cases} 0 & p \mid a\\ 1 & a是p的二次剩余\\ -1 & a不是p的二次剩余 \end{cases}$.
${p-1 \over (2, p-1)}={p-1 \over 2}$.
$1^2, 2^2, 3^2, \ldots, ({p-1 \over 2})^2$.
二次剩余数量$=$非二次剩余数量$={p-1 \over 2}$。
Euler判别法:$n^{p-1 \over 2} \equiv ({n \over p}) \pmod p, (n, p)=1$.
Gauss Lama:$n, 2n, 3n, \ldots, {p-1 \over 2}n \pmod p$的最小正剩余中大于${p-1 \over 2}$的个数为$m$,则$({n \over p})=(-1)^m$。
$({2 \over p})=(-1)^{p^2-1 \over 8}=(-1)^{[{p \over 2}]-[{p \over 4}]}$.
二次互反:$({q \over p})({p \over q})=(-1)^{{p-1 \over 2}{q-1 \over 2}}$.


Summary

$n=(a_ka_{k-1}a_{k-2}\ldots a_0)_p, m=(b_kb_{k-1}b_{k-2}\ldots b_0)_p, {n \choose m} \equiv \prod_{i=0}^k {a_i \choose b_i} \pmod p$.

$x^{\varphi(n)} \equiv 1 \pmod n, (x, n)=1$.
$x^n \equiv x^{n-\varphi(n)} \pmod n$.
$n$为素数,$x^{p-1} \equiv 1 \pmod p, p \not \mid x$.
$x^p \equiv x \pmod p$.

$x^k \equiv n \pmod p$无解或有$(k, p-1)$个解。
$x^k \equiv 1 \pmod 1$有$d$个解,$d=(k, p-1)=ak+b(p-1)$。
$x \Rightarrow p$的缩系,$x^k$有${p-1 \over (k, p-1)}$个不同的取值。
$\bmod p$的原根有$\varphi(p-1)$个。

$\lim \sum_{p=1}^{\infty} (1-{1 \over p})=0$.

$f(x) \equiv 0 \pmod {p^{\alpha}}$的解数$=f(x) \equiv 0 \pmod p$的解数$\Leftarrow \begin{cases} f(x) \equiv 0\\ f'(x) \equiv 0 \end{cases} \pmod p$无解。

$\begin{cases} b=2mn\\ a=n^2-m^2\\ c=n^2+m^2 \end{cases}$.

拉格朗日差值公式$f(x)=\sum_{i=1}^n \prod_{j \neq i} {(x-x_j) \over (x_i-x_j)}f(x_i)$.
牛顿差值公式$f(x)=a_0+a_1(x-x_1)+a_2(x-x_1)(x-x_2)+ \cdots +a_n \prod_{i=1}^n (x-x_i)$.

$f(x)=c_nx^n+c_{n-1}x^{n-1}+ \cdots +c_1x+c_0$可约$\Leftarrow p \mid c_0, p \mid c_1, \ldots, p \mid c_{n-1}, p \not \mid c_n, p^2 \not \mid c_0$.
$f(x)=c_nx^n+c_{n-1}x^{n-1}+ \cdots +c_1x+c_0$可约$\Leftarrow p \mid (c_0, c_1, \ldots, c_{n-1}), p \not \mid c_n, p^2 \not \mid c_0$.

$p^{\alpha} \mid \mid a-b, (a, p)=1, (b, p)=1, p \ge 3$.
$V_p(a^m-b^m)=\alpha+V_p(m)$.