BZOJ2142 礼物-扩展lucas

传送门

题意:

给定n及m个数a[1…m],再给定一个模数P,求$\sum_{i=1}^nC_{n-\sum_{j=1}^{i-1}a_j}^{a_i}$对P取模

($1≤n≤10_9,1≤m≤5$,设$ P=p_1^{c_1}\times p_2^{c_2}\times p_3^{c_3}\times …\times p_t^{c_t} $,$p_i$为质数,$1≤p_i^{c_i}≤10^5$)

Solution:

扩展lucas定理模板题

这个算法理解起来比较方便,但是写的时候就不是那么好写了…

我们把p拆成$p_1^{c_1}\times p_2^{c_2}\times p_3^{c_3}\times … \times p_t^{c_t}$的形式,这样我们只需要求出每个$C(n,m)%p_i^{c_i}$,再用CRT合并即可

而C(n,m)可以转化成阶乘及阶乘的逆元的乘积,由于扩欧求逆元要求互质,所以说我们计算阶乘时需要提取出和$p_i^{c_i}$互质的数

举一个其他博客用烂的栗子:

$p_i=3,c_i=2,n=19$

$n!=1\times 2\times 3\times 4\times 5\times 6\times 7\times 8\times ……\times 19$

$=[1\times 2\times 4\times 5\times 7\times 8\times … \times 16\times 17\times 19]\times (3\times 6\times 9\times 12\times 15\times 18)$

$=[1\times 2\times 4\times 5\times 7\times 8\times … \times 16\times 17\times 19]\times 3^6(1\times 2\times 3\times 4\times 5\times 6)$

发现后面的是$(\frac{n}{p_i})!$,于是递归即可。前半部分是以$p_i^{c_i}$为周期的$[1\times 2\times 4\times 5\times 7\times 8]=[10\times 11\times 13\times 14\times 16\times 17](\mod 9)$

最后会有孤立出来的19,可以证明孤立出的长度不超过$p_i^{c_i}$,暴力算即可

至于剩下的$3^6$之类的数,我们只要计算出n!,m!,(n-m)!里含有多少个$p_i$(不妨设a,b,c),那么a-b-c就是C(n,m)中p的个数,直接算一下就行。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include<cstdio>
#include<iostream>
using namespace std;
long long mod;
long long ans,sum;
long long a[10];
int n,m;
int nump,p;
long long tot;
int fast_pow(int a,int x,int mo)
{
int ans=1;
for (;x;x>>=1,a=1ll*a*a%mo)
if (x&1) ans=1ll*ans*a%mo;
return ans;
}
int mi(int n,int pi,int pk)
{
if (n==0) return 1;
int ans=1;
for (int i=2;i<=pk;i++)
if (i%pi) ans=1ll*ans*i%pk;
ans=fast_pow(ans,n/pk,pk);
for (int i=2;i<=n%pk;i++)
if (i%pi) ans=1ll*ans*i%pk;
return 1ll*ans*mi(n/pi,pi,pk)%pk;
}
void exgcd(int a,int b,long long &x,long long &y)
{
if (!b) x=1,y=0;
else exgcd(b,a%b,y,x),y-=a/b*x;
}
int inv(int a,int mo)
{
if (a==0) return 0;
long long x=0,y=0;
exgcd(a,mo,x,y);
x=(x%mo+mo)%mo;
if (x==0) x+=mo;
return x;
}
long long C(int n,int m,int pi,int pk)
{
if (m>n) return 0;
int a=mi(n,pi,pk);
int b=mi(m,pi,pk);
int c=mi(n-m,pi,pk);
int num=0;
for (int i=n;i;i/=pi) num+=i/pi;
for (int i=m;i;i/=pi) num-=i/pi;
for (int i=n-m;i;i/=pi) num-=i/pi;
int ans=1ll*a*inv(b,pk)%pk*inv(c,pk)%pk*fast_pow(pi,num,pk)%pk;
return 1ll*ans*(mod/pk)%mod*inv(mod/pk,pk)%mod;
}
int main()
{
scanf("%lld",&mod);
scanf("%d%d",&n,&m);
for (int i=1;i<=m;i++)
scanf("%d",&a[i]),tot+=a[i];
if (n<tot) {printf("Impossible\n");return 0;}
ans=1;
for (int i=1;i<=m;i++)
{
n-=a[i-1];
long long P=mod;sum=0;
for (int j=2;j*j<=P;j++)
{
if (P%j==0)
{
p=1;
while (P%j==0) P/=j,p*=j;
sum=(sum+C(n,a[i],j,p))%mod;
}
}
if (P>1) sum=(sum+C(n,a[i],P,P))%mod;
ans=(ans*sum)%mod;
}
printf("%lld\n",ans);
}
作者

Fizzmy

发布于

2018-01-18

更新于

2021-08-06

许可协议

评论