Lucas定理

　　之前的数论全都塞在一篇里，都快没法看了...现在我决定把它拆开。如果有一些比较零散的就先放在那里。

　　lx老师画风突变，给我们布置了大量毒瘤数论。于是我...开始做。

　　首先说一下$Lucas$定理：

　　就是这么一个很奇妙的公式，至于为什么...抱歉不会证。

　　所以它有什么用呢？主要是用于$n,m$非常大，但是$p$还不算很大的时候快速处理组合数取模。拿到这个式子，我们对于前一部分递归求解，后一部分就可以直接算。直接算也是需要技巧的!

　　·$p$真的非常小：开一个二维数组，用递推预处理组合数；

　　·$p$大约在$10^5$：用组合数的通项公式：

　　预处理阶乘以及阶乘的逆元，对于一般的题目这样就足够了，附一份代码：　

　　1 void init()
　　2 {
　　3     f[0]=inv[0]=1;
　　4     for (R i=1;i<p;++i)
　　5     {
　　6         f[i]=f[i-1]*i%p;
　　7         inv[i]=quick_pow(f[i],p-2);
　　8     }
　　9 }

　　·很多题解说Lucas只能处理$p<=10^5$的数据，然而并不是!观察发现，导致耗费时间最多的其实就是预处理逆元，达到了惊人的$O(PlogP)$。我们知道对于$1-n$的逆元是可以线性递推的，那么阶乘会不会也可以呢？答案是肯定的，感谢lx老师的毒瘤题让我学会了新知识。

　　1 inv[p-1]=qui(p-1,p-2);
　　2 for (R i=p-1;i>=1;--i)
　　3     inv[i-1]=1LL*inv[i]*i%p;

　　至于为什么能这么做...并不是非常清楚。

　　在cry大佬的教导下我发现这个东西真是非常显然了，思考一下逆元的本质，就会发现：

　　其实就是这样的简单，然而我也确实想不到。

　　下面粘一个第二种预处理的全代码：

　　卢卡斯定理：https://www.luogu.org/problemnew/show/P3807

 1 // luogu-judger-enable-o2
 2 # include <cstdio>
 3 # include <iostream>
 4 # define ll long long
 5 
 6 using namespace std;
 7 
 8 int T;
 9 ll f[100005],inv[100005];
10 ll n,m,p;
11 
12 ll c(ll n,ll m)
13 {
14     if(n<m) return 0;  
15     return (f[n]*inv[m]%p)*inv[n-m]%p;
16 }
17 
18 ll lucas(ll n,ll m)
19 {
20     if(m==0)
21         return 1;
22     else
23         return (ll)lucas(n/p,m/p)*c(n%p,m%p)%p;
24 }
25 
26 ll quick_pow(ll x,ll c)
27 {
28     ll s=1;
29     while (c)
30     {
31         if(c&1) s=s*x%p;
32         x=x*x%p;
33         c=c>>1;
34     }
35     return s%p;
36 }
37 
38 void firs()
39 {
40     f[0]=inv[0]=1;
41     for (int i=1;i<p;++i)
42     {
43         f[i]=f[i-1]*i%p;
44         inv[i]=quick_pow(f[i],p-2);
45     }
46 }
47 
48 int main()
49 {
50     scanf("%d",&T);
51     while (T--)
52     {
53         scanf("%lld%lld%lld",&n,&m,&p);
54         firs();
55         printf("%lld
",lucas(n+m,m));
56     }
57     return 0;
58 }

Lucas

　　这时候我们发现一个小问题：如果$P$不是质数呢？可以用一个非常毒瘤的$exlucas$，据说省选都不考的东西，然而lx老师也当做课后题布置了下来。

　　扩展lucas是有一种弱化版的，就是说P虽然不是质数，但是每个质因子也只出现一次，比如说古代猪文那道题：

　　[SDOI2010]古代猪文：https://www.luogu.org/problemnew/show/P2480

　　现在做这种题我NOIP是不是要凉...学姐说我应该看看费马小定理的题目，然后一查就查到了这个题;

　　题意概述：其实这题没法概述，不如贴个图：

　　所以这道题就是求这个式子的值啦。$n,g<=1,000,000,000$

　　这题的幂实在是比较大，所以快速幂也无能为力了，但我们还有别的方法。一般题目的模数就是用来限制答案大小的，然而这里的模数非常有用。根据欧拉定理，如果$n$为质数，$a^{(n-1)}space mod space n=0$,因为这道题的模数是个质数，所以一定互质，可以用这个定理来化简式子：(如果$G$正好等于模数就不互质了，直接特判一下输出0)。

　　这样只要把上面那个东西求出来就可以使用快速幂算$G$的幂次方了。$n$虽然是挺大，但是开完根号就不算非常大了，于是我们可以试一试在根号时间内分解因子。如果一个数是$n$的因子且它大于根号$n$，那么用$n$除以它必然可以得到一个小于$n$的因子，所以枚举小于根号$n$的数再用$n$去除，就可以得到$n$的所有因子，还是要特判一下完全平方数的情况防止加重了。

　　现在就只有一个问题了:

　　组合数取模!好的，lucas...

　　然而事情并没有那么简单，理智分析发现这回的模数不是质数了啊...所以写一个循环看看它有哪些质数因子：

　　其实还是比较好的，毕竟没有重复的因子，如果有重复因子那就更是麻烦到了一个新境界。所以怎么求呢？

　　首先对于这每一个质数因子进行$lucas$，求出来四个解：

　　因为要同时满足这四个方程，所以就用中国剩余定理合并一下答案即可.

  1 // luogu-judger-enable-o2
  2 # include <cstdio>
  3 # include <iostream>
  4 # include <cmath>
  5 # define LL long long
  6 
  7 const int mod=999911659;
  8 const int a[5]={2,3,4679,35617};
  9 int n,g,ans;
 10 int f[5][35699],ni[5][35699];
 11 int b[5];
 12 int x;
 13 
 14 void exgcd(int a,int b,int &x,int &y)
 15 {
 16     if(b==0)
 17     {
 18         x=1;
 19         y=0;
 20         return ;
 21     }
 22     exgcd(b,a%b,y,x);
 23     y-=(LL)a/b*x;
 24 }
 25 
 26 int qui_pow(int x,int n,int p) //求x^n%p 
 27 {
 28     int s=1;
 29     x=x%p;
 30     while (n)
 31     {
 32         if(n&1) s=(LL)x*s%p;
 33         x=(LL)x*x%p;
 34         n=n>>1;
 35     }
 36     return s%p;
 37 }
 38 
 39 void firs(int x)
 40 {
 41     f[x][0]=1;
 42     ni[x][0]=1;
 43     for (int i=1;i<a[x];++i)
 44     {
 45         f[x][i]=f[x][i-1]*i%a[x];
 46         ni[x][i]=qui_pow(f[x][i],a[x]-2,a[x]);
 47     }
 48 }
 49 
 50 int c(int n,int m,int x)
 51 {
 52     if(n<m) return 0;
 53     return ((LL)f[x][n]*ni[x][m]%a[x])*(LL)ni[x][n-m]%a[x];
 54 }
 55 
 56 int lucas(int n,int m,int x)
 57 {
 58     if(m==0)
 59         return 1;
 60     else
 61         return (LL)lucas(n/a[x],m/a[x],x)%a[x]*(LL)c(n%a[x],m%a[x],x)%a[x];
 62 }
 63 
 64 int crt()
 65 {
 66     int m,M=999911658,ans=0;
 67     for (int i=0;i<4;++i)
 68     {
 69         m=M/a[i];
 70         m%=a[i];
 71         int x,y;
 72         exgcd(m,a[i],x,y);
 73         x=(x%a[i]+a[i])%a[i];
 74         m=M/a[i];
 75         x=(LL)x*m%M;
 76         ans=(ans+(LL)x*b[i]%M)%M;
 77     }
 78     return ans;
 79 }
 80 
 81 signed main()
 82 {
 83     scanf("%d%d",&n,&g);
 84     if(mod==g)
 85     {
 86         printf("0");
 87         return 0;
 88     }
 89     for (int i=0;i<4;++i)
 90         firs(i);
 91     g=g%999911658;
 92     for (int i=1;i*i<=n;++i)
 93     {
 94         if (n%i) continue;
 95         for (int j=0;j<4;++j)
 96             b[j]=lucas(n,i,j);
 97         ans=(crt()+ans)%999911658;
 98         if(i*i!=n)
 99         {
100             for (int j=0;j<4;++j)
101                 b[j]=lucas(n,n/i,j);
102             ans=(crt()+ans)%999911658;
103         }
104     }
105     printf("%d",qui_pow(g,ans,999911659));
106     return 0;
107 }

古代猪文

　　礼物:https://www.luogu.org/problemnew/show/P2183

　　题意概述:

　　Exlucas板子题,几乎就是一路照着题解改过来的,等联赛完了再复习吧.

  1 # include <cstdio>
  2 # include <iostream>
  3 # define R register int
  4 # define ll long long
  5 
  6 using namespace std;
  7 
  8 ll n,m,p,sum;
  9 int h;
 10 ll w[7],a[12];
 11 ll pr[12],po[12];
 12 ll ans=1;
 13 
 14 ll qui (ll a,ll b,ll p)
 15 {
 16     ll s=1;
 17     while (b)
 18     {
 19         if(b&1LL) s=s*a%p;
 20         a=a*a%p;
 21         b=b>>1LL;
 22     }
 23     return s;
 24 }
 25 
 26 void ex_gcd (ll a,ll b,ll &x,ll &y)
 27 {
 28     if(!b) { x=1LL,y=0LL; return ; }
 29     ex_gcd(b,a%b,y,x);
 30     y-=a/b*x;
 31 }
 32 
 33 ll inv (ll a,ll p)
 34 {
 35     ll x,y;
 36     ex_gcd(a,p,x,y);
 37     return (x%p+p)%p;
 38 }
 39 
 40 ll crt ()
 41 {
 42     ll ans=0;
 43     ll x,y,m,pri;
 44     for (R i=1;i<=h;++i)
 45     {
 46         pri=po[i];
 47         m=p/pri;
 48         m%=pri;
 49         ex_gcd(m,pri,x,y);
 50         x=(x%pri+pri)%pri;
 51         m=p/pri;
 52         x=x*m%p;
 53         ans=(ans+x*a[i])%p;
 54     }
 55     return ans;
 56 }
 57 
 58 ll fac (ll n,ll Pr,ll Po)
 59 {
 60     if(n==0) return 1;
 61     ll ans=1;
 62     for (R i=2;i<=Po;++i)
 63         if(i%Pr) ans=ans*i%Po;
 64     ans=qui(ans,n/Po,Po);
 65     for (R i=2;i<=n%Po;++i)
 66         if(i%Pr) ans=ans*i%Po;
 67     return ans*fac(n/Pr,Pr,Po)%Po;
 68 }
 69 
 70 ll C (ll n,ll m,ll Pr,ll Po)
 71 {
 72     if(m>n) return 0;
 73     ll t=0;
 74     ll a=fac(n,Pr,Po),b=fac(m,Pr,Po),c=fac(n-m,Pr,Po);
 75     for (ll i=n;i;i/=Pr) t+=i/Pr;
 76     for (ll i=m;i;i/=Pr) t-=i/Pr;
 77     for (ll i=n-m;i;i/=Pr) t-=i/Pr;
 78     return qui(Pr,t,Po)*a*inv(b,Po)%Po*inv(c,Po)%Po;
 79 }
 80 
 81 ll Lucas (ll n,ll m)
 82 {
 83     for (R i=1;i<=h;++i)
 84         a[i]=C(n,m,pr[i],po[i])%p;
 85     return crt();
 86 }
 87 
 88 void div (ll p)
 89 {
 90     for (R i=2;i*i<=p;++i)
 91     {
 92         if(p%i==0)
 93         {
 94             ++h;
 95             pr[h]=i;
 96             po[h]=1;
 97             while (p%i==0) po[h]*=i,p/=i;
 98         }
 99     }
100     if(p>1) pr[++h]=p,po[h]=p;
101 }
102 
103 int main()
104 {
105     scanf("%lld",&p),div(p);
106     scanf("%lld%lld",&n,&m);
107     for (R i=1;i<=m;++i)
108         scanf("%lld",&w[i]),sum+=w[i];
109     if(sum>n) { puts("Impossible"); return 0; }
110     for (R i=1;i<=m;++i)
111         ans=(ans*Lucas(n,w[i]))%p,n-=w[i];
112     printf("%lld",ans);
113     return 0;
114 }

礼物

　　combination：https://www.lydsy.com/JudgeOnline/problem.php?id=2982

　　题意概述：

　　小清新题目，纯属模板；即使这样我也WA了两次，现在真是手残到一定境界了。

 1 # include <cstdio>
 2 # include <iostream>
 3 # define R register int
 4 
 5 using namespace std;
 6 
 7 const int p=10007;
 8 const int maxn=10010;
 9 int T;
10 int n,m;
11 int f[maxn],inv[maxn];
12 
13 int qui (int a,int b)
14 {
15     int s=1;
16     while (b)
17     {
18         if(b&1) s=s*a%p;
19         a=a*a%p;
20         b>>=1;
21     }
22     return s;
23 }
24 
25 void init()
26 {
27     f[0]=inv[0]=1;
28     for (R i=1;i<=p;++i)
29     {
30         f[i]=f[i-1]*i%p;
31         inv[i]=qui(f[i],p-2);
32     }
33 }
34 
35 int c (int n,int m)
36 {
37     if(n<m) return 0;
38     return (f[n]*inv[m]%p)*inv[n-m]%p;
39 }
40 
41 int Lucas (int n,int m)
42 {
43     if(m==0) return 1;
44     return Lucas(n/p,m/p)*c(n%p,m%p)%p;
45 }
46 
47 int main()
48 {
49     scanf("%d",&T);
50     init();
51     while (T--)
52     {
53         scanf("%d%d",&n,&m);
54         printf("%d
",Lucas(n,m));
55     }
56     return 0;
57 }

combination

　　超能粒子炮·改：https://www.lydsy.com/JudgeOnline/problem.php?id=4591

　　题意概述：

　　$n,k<=10^{18}$

　　首先根据$Lucas$打一个表.显然$np,n\%p$是不会变的,所以我们只看$i$的部分;

　　利用和式求和的方法,把相同的归类到一起,就可以迭代求解了:

 1 # include <cstdio>
 2 # include <iostream>
 3 # include <cstring>
 4 # define ll long long
 5 # define R register int
 6 
 7 using namespace std;
 8 
 9 const int p=2333;
10 int T;
11 ll n,k,t,ans,x;
12 int f[3000],inv[3000];
13 int C[p+1][p+1];
14 
15 ll qui (ll a,ll b)
16 {
17     ll s=1;
18     while (b)
19     {
20         if(b&1) s=s*a%p;
21         a=a*a%p;
22         b=b>>1;
23     }
24     return s%p;
25 }
26 
27 void init()
28 {
29     f[0]=inv[0]=1;
30     for (int i=1;i<=p;++i)
31     {
32         f[i]=f[i-1]*i%p;
33         inv[i]=qui(f[i],p-2);
34     }
35     C[0][0]=1;
36     for (R i=1;i<=p;++i)
37     {
38         C[i][0]=1;
39         for (R j=1;j<=i;++j)
40             C[i][j]=(C[i-1][j]+C[i-1][j-1])%p;
41     }
42     for (R i=0;i<=p;++i)
43         for (R j=1;j<=p;++j)
44             C[i][j]=(C[i][j-1]+C[i][j])%p;
45 }
46 
47 ll c (ll n,ll m)
48 {
49     if(n<m) return 0;
50     return (f[n]*inv[m]%p)*inv[n-m]%p;
51 }
52 
53 ll luc (ll n,ll m)
54 {
55     if(m==0) return 1;
56     else      return (ll)luc(n/p,m/p)*c(n%p,m%p)%p;
57 }
58 
59 int Lucas (ll n,ll k)
60 {
61     if(k<0) return 0;
62     if(n==0||k==0) return 1;
63     if(n<=p&&k<=p)
64         return C[n][k]; 
65     return ((ll)C[n%p][p-1]*Lucas(n/p,k/p-1)%p+luc(n/p,k/p)*C[n%p][k%p]%p)%p;
66 }
67 
68 int main()
69 {
70     scanf("%d",&T);
71     init();
72     while (T--)
73     {
74         scanf("%lld%lld",&n,&k);
75         k=min(k,n);
76         printf("%d
",Lucas(n,k));
77     }
78     return 0;
79 }

超能粒子炮·改

---shzr