lucas 以及 扩展

lucas定理

  1. 广义表达:$C_{m}^n \equiv C_{m_i}^{n_i} \pmod{p}$

    其中$m=m_kp^k+m_{k-1}p^{k-1}+m_{k-2}p^{k-2}+\cdots +m_1p+m_0$

    $n=n_kp^k+n_{k-1}p^{k-1}+n_{k-2}p^{k-2}+\cdots +n_1p+n_0$


  2. 实质上用它来对大的组合数取模,即$C_m^{n} \equiv C_{m/p}^{n/p} \times C_{m\%p}^{n\%p} \pmod p$

  3. 证明:

    首先$C_p^{n}=\frac{p(p-1)(p-2)\cdots(p-n+1)}{n!}$

    当$p$为质数的时候,可见$C_p^n \mod p=0$(除了$n=0$和$n=p$的情况)

    对于$(1+x)^p=\sum_{r=0}^p C_p^r x^r$,可知对于$r \in [1,p-1]$时该项都模$p$为$0$

    故$(1+x)^p \equiv 1+x^p \pmod p$

    进而有$(1+x)^{p^i} \equiv 1+x^{p^i} \pmod p$


    将$m$写成$p$进制的形式,则$m=\sum_{i=0}^k m_ip^i$,其中$m_i \in [0,p-1]$

    基于上述结论有

    $$\sum _{n=0}^m C_m^n x^n =(1+x)^m =\prod_{i=0}^k ((1+x)^{p^i})^{m_i} \equiv \prod _{i=0}^k (1+x^{p^i})^{m_i} =\prod_{i=0}^k (\sum_{n_i=0}^{m_i}C_{m_i}^{n_i}x^{n_ip^i}) $$

    $C_n^m$中$m>n$的情况为0,而$$\prod_{i=0}^k (\sum_{n_i=0}^{m_i}C_{m_i}^{n_i}x^{n_ip^i})=\prod_{i=0}^k (\sum_{n_i=0}^{m_i}C_{m_i}^{p-1}x^{n_ip^i})=\sum _{n=0}^{m}(\prod_{i=0}^k C_{m_i}^{n_i})x^n$$

    所以$$ \sum_{n=0}^m C_m^n x^n =\sum _{n=0}^{m}(\prod_{i=0}^k C_{m_i}^{n_i})x^n$$

    即$C_m^n=\prod_{i=0}^k C_{m_i}^{n_i}$,得证。

  4. 代码

    1
    2
    3
    4
    5
    6
    7
    8
    int Lucas(int n,int m)
    {
    if(n<m) return 0;
    else {
    if(n<=mo && m<=mo) return getC(n,m);
    else return 1ll*Lucas(n/mo,m/mo)*getC(n%mo,m%mo)%mo;
    }
    }

    其中$getC$可以是求小组合数的函数,也可以是直接预处理好的组合数(在模数较小的情况下)。

对lucas的扩展

  1. 问题:求$C_n^m \mod P$其中$P$为合数。

  2. 求解:设答案为$ans$

    则有

    $$ \begin{cases}ans \equiv c_1 \pmod {p_1^{k_1}}\\ ans \equiv c_2 \pmod {p_2^{k_1}} \\ \cdots \\ ans \equiv c_i \pmod {p_i^{k_i}} \end{cases}$$

    其中$c_i=C_n^m \mod p_i^{k_i}$

    由于是不同的质因数,所以可以求出对各个余数的解以后用中国剩余定理进行合并。

    求解一个组合数$C_n^m$对$p_i^{k_i}$的过程可以分别求分子和分母,然后再进行运算,即,求出$n!\mod p_i^{k_i}$,$m!\mod p_i^{k_i}$,$(n-m)!\mod p_i^{k_i}$后合并。注意分子可以含有大于$k_i$个$p_i$,这样的话求出来就直接为0了,而除掉分母的时候是可能会除掉一些$p_i$的;同时分母本身也可能有含有$p_i$的约数,于是统一把$p_i$移除来,不乘进答案,接下来用分子有的$p_i$的质数减掉分母的,然后乘上去即可。不可能出现减掉后为负数的情况,因为那样意味着组合数除不尽;有可能出现减掉后仍然大于$k_i$的情况,那么这个组合数在模$p_i^{k_i}$下为$0$。

    具体地,对于求$n! \mod p_i^{k_i}$可以分为几个部分求解。将每个$p_i$的倍数提出一个$p_i$来,组成$p_i^{\lfloor \frac{n}{p_i} \rfloor}$;对于剩余的部分,容易发现除掉了$p_i$的数(原来是$p_i$的倍数)组成了一个新的阶乘,即$\lfloor \frac{n}{p_i} \rfloor !$,这个我们可以递归求解。剩下的数形如$1,2,\cdots,p-1,p+1,p+2,\cdots,2p-1 \cdots$(因为$p_i$的倍数已经通过前两步瓜分到阶乘和次幂里去了),由于$a \times b \mod p=(a \mod p) \times (b \mod p)$,易得$\prod _{i=1}^{p_i^{k_i}-1}i[i \mod p \neq 0]$与$\prod _{i=p_i^{k_i}+1}^{2p_i^{k_i}-1}i[i \mod p \neq 0] $在$\mod p_i^{k_i}$下是相等的,那么只需求一个$p^k$的长度的乘积在模$p_i^{k_i}$下的结果即可,而这个阶乘里包含了$ \lfloor \frac{n}{p_i^{k_i}} \rfloor$个这个余数,快速幂即可;还有剩下来的、没有达到$p_i^{k_i}$的一段数,由于长度不会超过$p_i^{k_i}$,直接计算即可。

    求完所有的$c$就可以直接通过中国剩余定理来求$ans$了。

题目

poj 2891 Strange Way to Express Integers

扩展$crt$的板子…直接合并即可。

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
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
using namespace std;
const int maxk=1005;
typedef long long ll;
ll gcd(ll a,ll b) {return b==0?a:gcd(b,a%b);}
ll extgcd(ll a,ll b,ll &x,ll &y)
{
if(!b) {
x=1ll,y=0ll;return a;
}else {
int d=extgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
}
ll inv(ll tmp,ll p)
{
ll x,y;
extgcd(tmp,p,x,y);
x=(x+p)%p;
if(x==0) x+=p;
return x;
}
ll a[maxk],p[maxk];
int k;
int main()
{
while(scanf("%d",&k)==1) {
for(int i=1;i<=k;i++) {
scanf("%lld%lld",&p[i],&a[i]);
}
ll tmpm,tmpc,m1,m2,c1,c2;
bool can=true;
for(int i=2;i<=k;i++) {
c1=a[i-1],c2=a[i],m1=p[i-1],m2=p[i];
ll g=gcd(m1,m2);
if((c2-c1)%g!=0) {can=false;break;}
tmpm=m1*m2/g;
tmpc=(inv(m1/g,m2/g)*((c2-c1)/g))%(m2/g)*m1+c1;
tmpc=(tmpc+tmpm)%tmpm;
a[i]=tmpc,p[i]=tmpm;
}
if(!can) {
printf("-1\n");
}else printf("%lld\n",a[k]);
}
return 0;
}

bzoj 4591 [Shoi2015]超能粒子炮·改

题意:求$\sum_{i=0}^k C_{n}^i \mod 2333$

$$\begin{array}{lr} \sum_{i=0}^k C_{n}^i \\ = C_{n/p}^0 C_{n\%p}^0 + C_{n/p}^0 C_{n\%p}^1 +C_{n/p}^0 C_{n\%p}^2 + \cdots +C_{n/p}^0 C_{n\%p}^{2332} + C_{n/p}^1 C_{n\%p}^0 + C_{n/p}^1 C_{n\%p}^1 +C_{n/p}^1 C_{n\%p}^2 + \cdots +C_{n/p}^1 C_{n\%p}^{2332} + \cdots + \\ C_{n/p}^{\lfloor k/p \rfloor} C_{n\%p}^0 +\cdots + C_{n/p}^{\lfloor k/p \rfloor} C_{n\%p}^{k\%p} \\ = \sum_{i=0}^{\lfloor k/p \rfloor -1} (C_{n/p}^i \sum_{j=0}^{2332} C_{n\%p}^j) + C_{n/p}^{\lfloor k/p \rfloor} \times \sum_{i=0}^{k\%p} C_{n\%p}^i \\ =\sum_{j=0}^{2332} C_{n\%p}^j (\sum_{i=0}^{\lfloor k/p \rfloor -1} C_{n/p}^i ) + C_{n/p}^{\lfloor k/p \rfloor} \times \sum_{i=0}^{k\%p} C_{n\%p}^i\end{array}$$

对$\sum_{i=0}^{\lfloor k/p \rfloor -1} C_{n/p}^i $递归计算,对$C_{n/p}^{\lfloor k/p \rfloor}$直接lucas计算,对于小于$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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mo=2333;
int T;
ll N,K;
int c[mo+5][mo+5],sum[mo+5][mo+5];
void ini()
{
c[0][0]=1;
c[1][0]=c[1][1]=1;
for(int i=2;i<=mo;i++) {
c[i][0]=c[i][i]=1;
for(int j=1;j<i;j++) c[i][j]=(c[i-1][j-1]+c[i-1][j])%mo;
}
for(int i=0;i<=mo;i++) {
sum[i][0]=1;
for(int j=1;j<=mo;j++) sum[i][j]=(sum[i][j-1]+c[i][j])%mo;
}
}
int Lucas(ll n,ll m)
{
if(n<m) return 0;
else {
if(n<=mo && m<=mo) return c[n][m];
else return (1ll*Lucas(n/mo,m/mo)*c[n%mo][m%mo])%mo;
}
}
int cal(ll n,ll k)
{
if(n<=mo && k<=mo) return sum[n][k];
else {
return (1ll*sum[n%mo][2332]*cal(n/mo,k/mo-1)%mo+(1ll*Lucas(n/mo,k/mo)*sum[n%mo][k%mo])%mo)%mo;
}
}
int main()
{
ini();
scanf("%d",&T);
for(int z=0;z<T;z++) {
scanf("%lld%lld",&N,&K);
printf("%d\n",cal(N,K));
}
return 0;
}

bzoj 2142 礼物

比较裸的扩展,题解

bzoj 1951 [Sdoi2010]古代猪文

题意是求解$G^m \mod p$其中$m= \sum_{i=1}^{N} [N\%i==0]C_N^i$

这个就也比较裸,只是要套用一层费马小定理

根据质数循环节对指数只用求出模$\phi(P)-1$的结果就可以了,$P$是个质数,因而$P-1$为合数,故也考虑扩展$lucas$求解即可,然后快速幂搞一下。

需要注意的是当$G=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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include <bits/stdc++.h>
using namespace std;
const int P=999911659;
int ys[2000+5],al=0,N,G,c[5],rm[5],fac[4][36000+5];
map<int ,int > re;
const int pri[]={2,3,4679,35617};
void ysdiv(int x)
{
for(int i=1;i*i<=x;i++) {
if(x%i==0) {
ys[++al]=i;
if(x/i!=i) ys[++al]=x/i;
}
}
}
int quickpow(int x,int k,int mo)
{
x%=mo;
int ret=1;
while(k>0) {
if(k&1) ret=1ll*ret*x%mo;
x=1ll*x*x%mo;
k>>=1;
}
return ret;
}
inline int gcd(int a,int b) {return b==0?a:gcd(b,a%b);}
int extgcd(int a,int b,int &x,int &y)
{
if(b==0) {
x=1,y=0;
return a;
}else {
int d=extgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
}
int inv(int q,int p)
{
int x,y;
extgcd(q,p,x,y);
x=(x%p+p)%p;
if(!x) x+=p;
return x;
}
typedef pair<int ,int > pii;
#define fir first
#define sec second
#define MP make_pair
int Cnt(int x,int p)
{
if(x==0) return 0;
return Cnt(x/p,p)+x/p;
}
void ini()
{
for(int id=0;id<4;id++) {
fac[id][0]=1;
for(int j=1;j<pri[id];j++) fac[id][j]=1ll*fac[id][j-1]*j%pri[id];
}
}
pii Cal(int x,int id)
{
if(x==1 || x==0) return MP(1,0);
int p=pri[id];
int ret=1,del=Cnt(x,p);
ret=1ll*ret*quickpow(fac[id][p-1],x/p,p)%p;
pii nw=Cal(x/p,id);
ret=1ll*ret*nw.fir%p;
ret=1ll*ret*fac[id][x%p]%p;
return MP(ret,del);
}
int getC(int n,int m,int id)
{
if(n<m) return 0;
else {
pii up=Cal(n,id),dw1=Cal(m,id),dw2=Cal(n-m,id);
up.sec-=dw1.sec; up.sec-=dw2.sec;
assert(up.sec>=0);
if(up.sec>0) return 0;
else {
int ret=1;
ret=1ll*ret*up.fir%pri[id];
ret=1ll*ret*inv(dw1.fir,pri[id])%pri[id];
ret=1ll*ret*inv(dw2.fir,pri[id])%pri[id];
return ret;
}
}
}
int Lucas(int n,int m,int id)
{
if(n<m) return 0;
else {
if(n<pri[id] || m<pri[id]) return getC(n,m,id);
else return 1ll*Lucas(n/pri[id],m/pri[id],id)*1ll*getC(n%pri[id],m%pri[id],id)%pri[id];
}
}
int crt()
{
for(int i=0;i<4;i++) rm[i]=1;
for(int i=0;i<4;i++)
for(int j=0;j<4;j++) {
if(i==j) continue;
rm[j]=1ll*rm[j]*pri[i]%(P-1);
}
int ret=0;
for(int i=0;i<4;i++) {
ret=(ret+1ll*c[i]*rm[i]%(P-1)*1ll*inv(rm[i],pri[i])%(P-1))%(P-1);
}
return ret;
}
int main()
{
scanf("%d%d",&N,&G);
if(G==P) {printf("0\n");exit(0);}
ysdiv(N);
int M=0;
ini();
for(int i=1;i<=al;i++) {
for(int j=0;j<4;j++) c[j]=Lucas(N,ys[i],j);
M=(M+crt())%(P-1);
}
int res=1;
res=quickpow(G,M,P);
printf("%d\n",res);
return 0;
}