题目描述
组合数 $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