题目描述
组合数 \(C_n^m\) 表示的是从 \(n\) 个互不相同的物品中选出 \(m\) 个物品的方案数。举个例子， 从 \((1, 2, 3)\) 三个物品中选择两个物品可以有 \((1, 2)\)，\((1, 3)\)，\((2, 3)\) 这三种选择方法。根据组合数的定义，我们可以给出计算组合数 \(C_n^m\) 的一般公式：
\[C_n^m = \frac {n!} {m! \ (n - m)!} \]其中 \(n! = 1 \times 2 \times \cdots \times n\)。（特别地，当 \(n = 0\) 时，\(n! = 1\)；当 \(m > n\) 时，\(C_n^m = 0\)。）
小葱在 NOIP 的时候学习了 \(C_i^j\) 和 \(k\) 的倍数关系，现在他想更进一步，研究更多关于组合数的性质。小葱发现，\(C_i^j\) 是否是 \(k\) 的倍数，取决于 \(C_i^j \bmod k\) 是否等于 \(0\)，这个神奇的性质引发了小葱对 \(\mathrm{mod}\) 运算（取余数运算）的兴趣。现在小葱选择了是四个整数 \(n, p, k, r\)，他希望知道
\[\left( \sum_{i = 0}^\infty C_{nk}^{ik + r} \right) \bmod p, \]即
\[\left( C_{nk}^{r} + C_{nk}^{k + r} + C_{nk}^{2k + r} + \cdots + C_{nk}^{(n - 1)k + r} + C_{nk}^{nk + r} + \cdots \right) \bmod p \]的值。
输入格式
第一行有四个整数 \(n, p, k, r\)，所有整数含义见问题描述。
输出格式
一行一个整数代表答案。
数据范围与提示
对于 \(30\%\) 的测试点，\(1 \leq n, k \leq 30\)，\(p\) 是质数；
对于另外 \(5\%\) 的测试点，\(p = 2\)；
对于另外 \(5\%\) 的测试点，\(k = 1\)；
对于另外 \(10\%\) 的测试点，\(k = 2\)；
对于另外 \(15\%\) 的测试点，\(1 \leq n \leq 10^3, 1 \leq k \leq 50\)，\(p\) 是质数；
对于另外 \(15\%\) 的测试点，\(1 \leq n \times k \leq 10^6\)，\(p\) 是质数；
对于另外 \(10\%\) 的测试点，\(1 \leq n \leq 10^9, 1 \leq k \leq 50\)，\(p\) 是质数；
对于 \(100\%\) 的测试点，\(1 \leq n \leq 10^9, 0 \leq r < k \leq 50, 2 \leq p \leq 2^{30} - 1\)。
大概是我经验少了，完全看不出来是 DP 呢daze。
乍一看题目数据范围大的离谱： \(\infty\) ?!
好怪哦，再看一眼：
\[\left( C_{nk}^{r} + C_{nk}^{k + r} + C_{nk}^{2k + r} + \cdots + C_{nk}^{(n - 1)k + r} + C_{nk}^{nk + r} + \cdots \right) \bmod p \]哦，这就没事了daze。
从 \(C_{nk}^{nk + r} \bmod p\) 开始，剩下的项就都是 \(0\) 了。
要求的就是
\[\left( \sum_{i = 0}^{n - 1} C_{nk}^{ik + r} \right) \bmod p 。 \]于是就可以预处理阶乘和阶乘逆元，暴力求每一项组合数，（期望）可以拿到 \(n\) 比较小的暴力分。
对于大一点的组合数，可以用 Lucas 定理骗点分。
然而 \(p\) 有可能很大，所以还是 RE 很多点。
但还是能拿 \(60\) ，很多了（自我安慰）。
#include<iostream>
#include<cstdio>
#define int long long
using namespace std;
const int N=1e6+5;
int n,mod,k,r;
int fact[N],inv[N];
inline int qpow(int a,int p){
int ret=1;
while(p){
if(p&1)(ret*=a)%=mod;
(a*=a)%=mod,p>>=1;
}
return ret;
}
inline int C(int n,int m){
if(n<0||m<0||n<m)return 0;
if(n>=mod||m>=mod)return C(n/mod,m/mod)*C(n%mod,m%mod)%mod;//Lucas 定理
return fact[n]*inv[m]%mod*inv[n-m]%mod;
}
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
cin>>n>>mod>>k>>r;
fact[0]=1;
for(int i=1;i<min(mod,N);i++)fact[i]=fact[i-1]*i%mod;
inv[min(mod,N)-1]=qpow(fact[min(mod,N)-1],mod-2);
for(int i=min(mod,N)-2;~i;i--)inv[i]=inv[i+1]*(i+1)%mod;//预处理阶乘和阶乘逆元
int ans=0;
for(int i=r;i<=n*k;i+=k)(ans+=C(n*k,i))%=mod;//暴力求答案
cout<<ans<<endl;
return 0;
}
时间复杂度是 \(\Theta(魔理沙)\) 。
暴力打完思考一下正解。
\(n\) 的范围也太大了点吧，如果要求每一项，至少也要 \(\Theta(n)\) ，毫无前途。
那我们试试推下式子？
……
以我目前的水平显然推不了。
但是 \(k\) 和 \(r\) 的范围好小啊，只有 \(50\) ，而且 \(0 \leq r < k\) ？什么意思？有必要吗？
好像想到了什么……
\(r\) 像是余数？
\(k\) 像是除数？
好像的确是这样，我们换个角度，考虑这个东西的组合意义，
要求的东西就是在 \(nk\) 个球中选 \(j\) 个，且 \(j \bmod k = r\)的方案数。
那我们试试 DP ，设 \(f[i][j]\) 表示考虑前 i 个球，选的球数 \(m\) 满足 \(m \bmod k = j\) 的方案数。
想一想怎么转移，每个球有选或不选两种选择。
选的话就是 \(f[i - 1][( j - 1 + k ) \% k]\) ，不选就是 \(f[i - 1][j]\) 。
所以 \(f[i][j] = f[i - 1][( j - 1 + k ) \% k] + f[i - 1][j]\) 。
直接转移会 MLE ，但每次转移只与 \(f[i - 1]\) 有关，显然可以滚动。
#include<iostream>
#include<cstdio>
#define int long long
using namespace std;
const int K=55;
int n,mod,k,r;
int f[2][K];
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
cin>>n>>mod>>k>>r;
f[0][0]=1;
int cur=0;
for(int i=1;i<=n*k;i++){
cur^=1;
for(int j=0;j<k;j++)
f[cur][j]=(f[cur^1][j]+f[cur^1][(j-1+k)%k])%mod;
}
cout<<f[cur][r]<<endl;
return 0;
}
时间复杂度是 \(\Theta(nk^2)\) 的。
依然只能拿 \(60\) \kk。
虽然但是，这份代码相当有前途。
重新审视一下数据范围：
\(n\) 大，\(k\) 小。
\(n\) 很大，\(k\) 很小。
\(n\) 非常大，\(k\) 非常小。
对于 DP 来说，\(10^9\) 的数据范围，要么数位DP，要么倍增或者矩阵快速幂。
对于 DP 来说，\(50\) 的数据范围，要么状压DP，要么矩阵快速幂。
显然，可以考虑矩阵快速幂。
矩阵快速幂的本质是优化转移过程，非常适合优化像这样每次转移只与 \(i - 1\) 有关并且只有加法的状态转移方程。
设计转移矩阵 \(C\) 时，应该让状态矩阵 \(A\) 满足 \(A_{i-1} \times C = A_i\) 。
所以 \(C\) 应该为
\[\left[ \begin{matrix} 1 & 0 & 0 & \ldots & 0 & 0 & 1\\ 1 & 1 & 0 & \ldots & 0 & 0 & 0\\ 0 & 1 & 1 & \ldots & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & \dots & 1 & 1 & 0\\ 0 & 0 & 0 & \dots & 0 & 1 & 1\\ \end{matrix} \right] \]但 \(k = 1\) 是边界会有些问题，特判一下就好。
然后直接矩乘就好了。
代码实现
#include<iostream>
#include<cstdio>
#include<cstring>
#define int long long
using namespace std;
const int K=55;
int n,mod,k,r;
int f[2][K];
struct matrix{
int n,m;
int a[K][K];
matrix(){memset(a,0,sizeof(a));n=0,m=0;}
friend inline matrix operator*(const matrix&a,const matrix&b){
matrix c;
c.n=a.n,c.m=b.m;
for(int i=1;i<=a.n;i++)
for(int j=1;j<=b.m;j++)
for(int k=1;k<=a.m;k++)
c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%mod)%mod;
return c;
}
};//普通的矩阵模板
inline matrix qpow(matrix a,int p){//普通的矩阵快速幂
matrix ret;
ret.n=a.n,ret.m=a.m;
for(int i=1;i<=a.n;i++)ret.a[i][i]=1;
while(p){
if(p&1)ret=ret*a;
a=a*a,p>>=1;
}
return ret;
}
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
cin>>n>>mod>>k>>r;
matrix C,I;
C.n=k,C.m=k;
if(k==1)C.a[1][1]=2;//特判
else{
C.a[1][1]=1,C.a[1][k]=1;
for(int i=2;i<=k;i++)
C.a[i][i-1]=1,C.a[i][i]=1;//构造转移矩阵
}
I.n=1,I.m=k;
I.a[1][1]=1;//初始矩阵
matrix B=I*qpow(C,n*k);//答案矩阵
cout<<B.a[1][r+1]<<endl;
return 0;
}
//iictiw-marisa