卢卡斯定理：
C
n
m
=
C
n
/
p
m
/
p
∗
C
n
%
p
m
%
p
(
m
o
d
p
)
C_n^m=C_{n/p}^{m/p}*C_{n\%p}^{m\%p}(mod\ p)
Cnm=Cn/pm/p∗Cn%pm%p(mod p)
其中
p
p
p是质数
证明没有去搜，想了解请自行搜索吧。我觉得式子不难记，所以当个结论记就好了。复杂度我也没仔细考虑过，但是据说可以做模数是
1
e
5
1e5
1e5级别的问题。
扩展卢卡斯定理
若
p
p
p不是质数怎么办呢？
这时候就需要扩展卢卡斯定理了。
我们对
p
p
p进行质因数分解，令
p
=
p
1
k
1
∗
p
2
k
2
.
.
.
p
n
k
n
p=p_1^{k_1}*p_2^{k_2}...p_n^{k_n}
p=p1k1∗p2k2...pnkn
那么我们如果能求出
C
n
m
C_n^m
Cnm对于每一个
p
i
k
i
p_i^{k_i}
piki的结果，由于每一个
p
i
k
i
p_i^{k_i}
piki都是同一质因子的若干次方，所以不同质因子的若干次方仍然是互质的，那么我们可以用中国剩余定理对方程组合并来得到最终的结果。
那么我们考虑如何求出
C
n
m
C_n^m
Cnm对于每一个
p
i
k
i
p_i^{k_i}
piki。
我们知道，
C
n
m
=
n
!
m
!
(
n
−
m
)
!
C_n^m=\frac{n!}{m!(n-m)!}
Cnm=m!(n−m)!n!，那么我们需要求出
n
!
%
p
i
k
i
n!\%p_i^{k_i}
n!%piki、
m
!
%
p
i
k
i
m!\%p_i^{k_i}
m!%piki和
(
n
−
m
)
!
%
p
i
k
i
(n-m)!\%p_i^{k_i}
(n−m)!%piki，然后再利用逆元即可。
对于
n
!
%
p
i
k
i
n!\%p_i^{k_i}
n!%piki，我们把
n
!
n!
n!拆开分解成三个部分，第一部分是所有不含因子
p
i
p_i
pi的数的乘积，我们发现这个是每
p
i
k
i
p_i^{k_i}
piki个数取模后的结果相同，也就是有循环节的，那么我们算出前
p
i
k
i
p_i^{k_i}
piki个数中没有被提出因子
p
i
p_i
pi的数在模
p
i
k
i
p_i^{k_i}
piki意义下的乘积记为
x
x
x，这个乘积要算
⌊
n
p
i
k
i
⌋
\lfloor\frac{n}{p_i^{k_i}}\rfloor
⌊pikin⌋次，那么之后用快速幂求出
x
⌊
n
p
i
k
i
⌋
x^{\lfloor\frac{n}{p_i^{k_i}}\rfloor}
x⌊pikin⌋，然后最后剩余的一点零散的部分暴力求解即可，因为零散的部分一点不超过
p
i
k
i
p_i^{k_i}
piki个数。第二部分是提出了
x
x
x个
p
i
p_i
pi，那么就要乘
p
i
x
p_i^{\ \ x}
pi x，可以快速幂求解。第三部分是提出来一个
p
i
p_i
pi后的那些数，而这些数又是一个新的阶乘，可以递归下去求解。
看到网上都举了那个
n
=
19
n=19
n=19，
p
=
3
p=3
p=3，
k
=
2
k=2
k=2的例子，我觉得还是很不错的，所以我也用一下。
对于求
19
!
%
3
2
19!\%3^2
19!%32，我们的操作过程如下：
n
!
=
1
∗
2
∗
3
∗
4
∗
5
∗
6
∗
7
∗
8
∗
9
∗
10
∗
11
∗
12
∗
13
∗
14
∗
15
∗
16
∗
17
∗
18
∗
19
n!=1∗2∗3∗4∗5∗6∗7∗8∗9∗10∗11∗12∗13∗14∗15∗16∗17∗18∗19
n!=1∗2∗3∗4∗5∗6∗7∗8∗9∗10∗11∗12∗13∗14∗15∗16∗17∗18∗19
=
(
1
∗
2
∗
∗
4
∗
5
∗
∗
7
∗
8
∗
∗
10
∗
11
∗
13
∗
14
∗
16
∗
17
∗
19
)
∗
(
3
∗
6
∗
9
∗
12
∗
15
∗
18
)
=(1∗2∗∗4∗5∗∗7∗8∗∗10∗11∗13∗14∗16∗17∗19)*(3*6*9*12*15*18)
=(1∗2∗∗4∗5∗∗7∗8∗∗10∗11∗13∗14∗16∗17∗19)∗(3∗6∗9∗12∗15∗18)
=
(
1
∗
2
∗
4
∗
5
∗
7
∗
8
∗
10
∗
11
∗
13
∗
14
∗
16
∗
17
∗
19
)
∗
3
6
∗
(
1
∗
2
∗
3
∗
4
∗
5
∗
6
)
=(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)∗3^6∗(1∗2∗3∗4∗5∗6)
=(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)∗36∗(1∗2∗3∗4∗5∗6)
实际上我们在计算时是看整个
C
n
m
C_n^m
Cnm中有多少个质因数
p
i
p_i
pi，计算
n
!
n!
n!中质因子
p
p
p的个数
x
x
x的公式为
x
=
⌊
n
p
⌋
+
⌊
n
p
2
⌋
+
⌊
n
p
3
⌋
+
.
.
.
x=\lfloor\frac{n}{p}\rfloor+\lfloor\frac{n}{p^2}\rfloor+\lfloor\frac{n}{p^3}\rfloor+...
x=⌊pn⌋+⌊p2n⌋+⌊p3n⌋+...
计算完之后再分别把
n
!
n!
n!，
m
!
m!
m!和
(
n
−
m
)
!
(n-m)!
(n−m)!的剩余部分求出来即可。
模板：bzoj3283运算器
注：
上面讲的是求
C
n
m
%
p
C_n^m\%p
Cnm%p的过程，代码写的是求
C
m
n
%
p
C_m^n\%p
Cmn%p的过程。代码是bzoj3283运算器的代码，其中第三个部分是扩展卢卡斯定理的部分，如果想学扩展BSGS，请点击这里。
代码有两种写法，区别在于CRT的合并上，第一种方法是用原来写CRT的写法来写，每一个式子都算好了再合并。另一种做法是因为我们最后所有的模数的乘积就是一开始给出的
p
p
p，所以我们可以一边求一边合并。
第一种代码：
#include <bits/stdc++.h>
using namespace std;
int t,opt;
long long p[50],a[50];//存储质因数p的t次方,a存储CRT要用的余数
map<long long,long long> mp;
long long ji;
inline long long ksm(long long x,long long y,long long mod)
{
long long res=1;
while(y)
{
if(y&1)
res=(res*x)%mod;
x=(x*x)%mod;
y>>=1;
}
return res;
}
inline void exgcd(long long a,long long b,long long &x,long long &y)
{
if(!b)
{
x=1;
y=0;
}
else
{
exgcd(b,a%b,y,x);
y-=a/b*x;
}
}
inline long long inv(long long a,long long b)
{
long long x=0,y=0;
exgcd(a,b,x,y);
x=(x%b+b)%b;
if(!x)
x+=b;
return x;
}
inline long long cal(long long n,long long x,long long mod)//递归计算除了x的若干次方之外的部分(第一、第三部分)
{
if(!n)
return 1;
long long ans=1;//提出来的那些不含因子x的乘积
if(n/mod)//有整块的
{
for(int i=1;i<=mod;++i)
{
if(i%x)//不含因子x
ans=ans*i%mod;
}
ans=ksm(ans,n/mod,mod);//有循环节，所以乘积用快速幂计算即可
}
for(int i=n/mod*mod+1;i<=n;++i)
{
if(i%x)
ans=ans*i%mod;
}
return ans*cal(n/x,x,mod)%mod;//当前的不含因子x的乘积乘以递归下去求的剩余阶乘部分的结果
}
//计算出对于每一个质数的若干次方取模后的结果
inline long long exlucas(long long m,long long n,long long x,long long mod)//x是当前质数
{
if(n>m)
return 0;
int cnt=0;
for(int i=m;i;i/=x)
cnt+=i/x;
for(int i=n;i;i/=x)
cnt-=i/x;
for(int i=m-n;i;i/=x)
cnt-=i/x;
long long ans=ksm(x,cnt,mod)*cal(m,x,mod)%mod*inv(cal(n,x,mod),mod)%mod*inv(cal(m-n,x,mod),mod)%mod;
return ans*(ji/mod)%ji*inv(ji/mod,mod)%ji;//CRT合并！
}
inline long long gcd(long long x,long long y)
{
return y?gcd(y,x%y):x;
}
int main()
{
scanf("%d",&t);
while(t--)
{
scanf("%d",&opt);
if(opt==1)
{
long long x,y,p;
scanf("%lld%lld%lld",&x,&y,&p);
printf("%lld\n",ksm(x,y,p));
}
else if(opt==2)
{
mp.clear();
long long aa,b,p;
int pd=0;
scanf("%lld%lld%lld",&aa,&b,&p);
if(b==1)
{
printf("0\n");
pd=1;
}
if(pd==1)
continue;
long long d=gcd(aa,p),t=1,k=0;
while(d!=1)
{
if(b%d)
{
printf("Math Error\n");
pd=1;
break;
}
++k;
b/=d;
p/=d;
t=(t*(aa/d))%p;
if(b==t)
{
printf("%lld\n",k);
pd=1;
break;
}
d=gcd(aa,p);
}
if(pd==1)
continue;
long long m=ceil(sqrt(p)),ans;
for(int j=0;j<=m;++j)
{
if(j==0)
{
ans=b%p;
mp[ans]=j;
continue;
}
ans=(ans*aa)%p;
mp[ans]=j;
}
long long x=ksm(aa,m,p);
ans=t;
for(int i=1;i<=m;++i)
{
ans=(ans*x)%p;
if(mp[ans])
{
x=i*m-mp[ans];
printf("%lld\n",x+k);
pd=1;
break;
}
}
if(!pd)
printf("Math Error\n");
}
else
{
long long x,y,mod;
scanf("%lld%lld%lld",&x,&y,&mod);
ji=mod;//mod就是CRT的那个因子乘积和！
if(mod==1)
{
printf("0\n");
continue;
}
memset(p,0,sizeof(p));
memset(a,0,sizeof(a));
int cnt=0;//记录质数个数
for(int i=2;i*i<=mod;++i)
{
if(mod%i==0)
{
p[++cnt]=1;
while(mod%i==0)
{
p[cnt]*=i;
mod/=i;//把当前因子除掉对看是否是后面数的倍数无影响
}
a[cnt]=exlucas(y,x,i,p[cnt]);
}
}
if(mod>1)
{
p[++cnt]=mod;
a[cnt]=exlucas(y,x,mod,mod);
}
long long res=0;
for(int i=1;i<=cnt;++i)
res=(res+a[i])%ji;
printf("%lld\n",res);
}
}
return 0;
}
第二种写法：
#include <bits/stdc++.h>
using namespace std;
int t,opt;
long long p[50],a[50];//存储质因数p的t次方,a存储CRT要用的余数
map<long long,long long> mp;
long long ji;
inline long long ksm(long long x,long long y,long long mod)
{
long long res=1;
while(y)
{
if(y&1)
res=(res*x)%mod;
x=(x*x)%mod;
y>>=1;
}
return res;
}
inline void exgcd(long long a,long long b,long long &x,long long &y)
{
if(!b)
{
x=1;
y=0;
}
else
{
exgcd(b,a%b,y,x);
y-=a/b*x;
}
}
inline long long inv(long long a,long long b)
{
long long x=0,y=0;
exgcd(a,b,x,y);
x=(x%b+b)%b;
if(!x)
x+=b;
return x;
}
inline long long cal(long long n,long long x,long long mod)//递归计算除了x的若干次方之外的部分
{
if(!n)
return 1;
long long ans=1;//提出来的那些不含因子x的乘积
if(n/mod)//有整块的
{
for(int i=1;i<=mod;++i)
{
if(i%x)//不含因子x
ans=ans*i%mod;
}
ans=ksm(ans,n/mod,mod);//有循环节，所以乘积用快速幂计算即可
}
for(int i=n/mod*mod+1;i<=n;++i)
{
if(i%x)
ans=ans*i%mod;
}
return ans*cal(n/x,x,mod)%mod;//当前的不含因子x的乘积乘以递归下去求的剩余阶乘部分的结果
}
//计算出对于每一个质数的若干次方取模后的结果
inline long long exlucas(long long m,long long n,long long x,long long mod)//x是当前质数
{
if(n>m)
return 0;
int cnt=0;
for(int i=m;i;i/=x)
cnt+=i/x;
for(int i=n;i;i/=x)
cnt-=i/x;
for(int i=m-n;i;i/=x)
cnt-=i/x;
long long ans=ksm(x,cnt,mod)*cal(m,x,mod)%mod*inv(cal(n,x,mod),mod)%mod*inv(cal(m-n,x,mod),mod)%mod;
return ans*(ji/mod)%ji*inv(ji/mod,mod)%ji;//CRT合并！
}
inline long long crt(int n)//CRT合并刚才得出的余数
{
long long m=1,ans=0;
for(int i=1;i<=n;++i)
m*=p[i];
for(int i=1;i<=n;++i)
{
long long x=0,y=0;
exgcd(m/p[i],p[i],x,y);
x=(x%p[i]+p[i])%p[i];
if(!x)
x+=p[i];
ans=(ans+a[i]*(m/p[i])*x%ji)%ji;
if(!ans)
ans+=ji;
}
return ans;
}
inline long long gcd(long long x,long long y)
{
return y?gcd(y,x%y):x;
}
int main()
{
scanf("%d",&t);
while(t--)
{
scanf("%d",&opt);
if(opt==1)
{
long long x,y,p;
scanf("%lld%lld%lld",&x,&y,&p);
printf("%lld\n",ksm(x,y,p));
}
else if(opt==2)
{
mp.clear();
long long aa,b,p;
int pd=0;
scanf("%lld%lld%lld",&aa,&b,&p);
if(b==1)
{
printf("0\n");
pd=1;
}
if(pd==1)
continue;
long long d=gcd(aa,p),t=1,k=0;
while(d!=1)
{
if(b%d)
{
printf("Math Error\n");
pd=1;
break;
}
++k;
b/=d;
p/=d;
t=(t*(aa/d))%p;
if(b==t)
{
printf("%lld\n",k);
pd=1;
break;
}
d=gcd(aa,p);
}
if(pd==1)
continue;
long long m=ceil(sqrt(p)),ans;
for(int j=0;j<=m;++j)
{
if(j==0)
{
ans=b%p;
mp[ans]=j;
continue;
}
ans=(ans*aa)%p;
mp[ans]=j;
}
long long x=ksm(aa,m,p);
ans=t;
for(int i=1;i<=m;++i)
{
ans=(ans*x)%p;
if(mp[ans])
{
x=i*m-mp[ans];
printf("%lld\n",x+k);
pd=1;
break;
}
}
if(!pd)
printf("Math Error\n");
}
else
{
long long x,y,mod;
scanf("%lld%lld%lld",&x,&y,&mod);
ji=mod;//mod就是CRT的那个因子乘积和！
if(mod==1)
{
printf("0\n");
continue;
}
memset(p,0,sizeof(p));
memset(a,0,sizeof(a));
int cnt=0;//记录质数个数
for(int i=2;i*i<=mod;++i)
{
if(mod%i==0)
{
p[++cnt]=1;
while(mod%i==0)
{
p[cnt]*=i;
mod/=i;//把当前因子除掉对看是否是后面数的倍数无影响
}
a[cnt]=exlucas(y,x,i,p[cnt]);
}
}
if(mod>1)
{
p[++cnt]=mod;
a[cnt]=exlucas(y,x,mod,mod);
}
printf("%lld\n",crt(cnt));
}
}
return 0;
}
6047
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
