题目传送门

求组合数的时候,如果模数p是质数,可以用卢卡斯定理解决。

但是卢卡斯定理仅仅适用于p是质数的情况。

当p不是质数的时候,我们就需要用扩展卢卡斯求解。

实际上,扩展卢卡斯=快速幂+快速乘+exgcd求逆元+质因数分解+crt合并答案+求阶乘,跟卢卡斯定理没什么关系......

如果把模数p分解成p1^k1*p2^k2*...*px^kx的形式,那么我们可以求出c(n,m)分别模每个pi^ki的结果,再用中国剩余定理合并即可。

每个pi^ki一定是互质的,所以用朴素crt就行。

根据组合数的定义,c(n,m)=(n!) / (m!*(n-m)!) ,所以我们只要能想办法求出阶乘,就能再利用exgcd求出逆元,进而求出组合数。

接下来唯一的问题就是怎么快速求出 x! 取模 pi^ki 的结果。

考虑如下的经典样例(据说来自popoqqq):(19!)%(3^2)

19!=1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19

先把其中的3的倍数提出来,因为求组合数的时候分子分母能约掉。

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)*(1*2*4*5*7*8)*(3*3*3*3*3*3)*(1*2*3*4*5*6)=(1*2*4*5*7*8)^2*19*(3^6)*(1*2*3*4*5*6)。

后面的6!部分可以递归求解,递归终点为0!=1。

3^6最后计算组合数的时候再处理。

那几个(1*2*4*5*7*8)显然是循环的,循环节长度小于pi^ki,可以暴力计算。

显然一共有(x/(pi^ki))个循环节,套个快速幂即可。

剩下的部分,即19,长度等于x%(pi^ki),也小于pi^ki,也可以暴力计算。

至此我们求出了阶乘。

求组合数的时候,考虑pi的倍数的影响。

分子分母分别计数相加减。

最后用crt合并即可。

 1 #include<cstdio>
 2 typedef long long ll;
 3
 4 ll n,m,p;
 5
 6 ll ksm(ll b,ll tp,ll mod)
 7 {
 8     ll ret=1;
 9     while(tp)
10     {
11         if(tp&1)ret=ret*b%mod;
12         b=b*b%mod;
13         tp>>=1;
14     }
15     return ret;
16 }
17
18 ll mul(ll a,ll b,ll mod)
19 {
20     ll ret=0;
21     while(b)
22     {
23         if(b&1)ret=(ret+a)%mod;
24         a=(a+a)%mod;
25         b>>=1;
26     }
27     return ret;
28 }
29
30 ll exgcd(ll a,ll b,ll &x,ll &y)
31 {
32     if(!b)
33     {
34         x=1;y=0;
35         return a;
36     }
37     ll t=exgcd(b,a%b,y,x);
38     y-=a/b*x;
39 }
40
41 ll inv(ll x,ll mod)
42 {
43     ll a,b;
44     exgcd(x,mod,a,b);
45     return (a%mod+mod)%mod;
46 }
47
48 ll fac(ll x,ll pi,ll pk)
49 {
50     if(!x)return 1;
51     ll ans=1;
52     for(ll i=2;i<=pk;i++)
53         if(i%pi)ans=ans*i%pk;
54     ans=ksm(ans,x/pk,pk);
55     for(ll i=2;i<=x%pk;i++)
56         if(i%pi)ans=ans*i%pk;
57     return ans*fac(x/pi,pi,pk)%pk;
58 }
59
60 ll c(ll cn,ll cm,ll pi,ll pk)
61 {
62     if(cm>cn)return 0;
63     ll up=fac(cn,pi,pk),d1=fac(cm,pi,pk),d2=fac(cn-cm,pi,pk);
64     ll cnt=0;
65     for(ll i=cn;i;i/=pi)cnt+=i/pi;
66     for(ll i=cm;i;i/=pi)cnt-=i/pi;
67     for(ll i=cn-cm;i;i/=pi)cnt-=i/pi;
68     return up*inv(d1,pk)%pk*inv(d2,pk)%pk*ksm(pi,cnt,pk)%pk;
69 }
70
71 ll crt(ll a,ll pk)
72 {
73     return a*inv(p/pk,pk)%p*(p/pk)%p;
74 }
75
76 int main()
77 {
78     scanf("%lld%lld%lld",&n,&m,&p);
79     ll tp=p,ans=0;
80     for(ll i=2;i*i<=p;i++)
81     {
82         if(tp%i)continue;
83         ll pk=1;
84         while(!(tp%i))tp/=i,pk*=i;
85         ans=(ans+crt(c(n,m,i,pk),pk))%p;
86     }
87     if(tp>1)ans=(ans+crt(c(n,m,tp,tp),tp))%p;
88     printf("%lld",(ans%p+p)%p);
89     return 0;
90 }

转载于:https://www.cnblogs.com/eternhope/p/9898494.html

[洛谷P4720] [模板] 扩展卢卡斯相关推荐

  1. 洛谷P4720 【模板】扩展卢卡斯

    P4720 [模板]扩展卢卡斯 题目背景 这是一道模板题. 题目描述 求 C(n,m)%P 其中 C 为组合数. 输入输出格式 输入格式: 一行三个整数 n,m,p ,含义由题所述. 输出格式: 一行 ...

  2. 洛谷 P4720 【模板】扩展卢卡斯定理/exLucas

    [模板]扩展卢卡斯定理/exLucas 题目背景 这是一道模板题. 题目描述 求 Cnmmodp{\mathrm{C}}_n^m \bmod{p}Cnm​modp 其中 C\mathrm{C}C 为组 ...

  3. 专题·树链剖分【including 洛谷·【模板】树链剖分

    初见安~~~终于学会了树剖~~~ [兴奋]当初机房的大佬在学树剖的时候我反复强调过:"学树剖没有前途的!!!" 恩.真香. 一.重链与重儿子 所谓树剖--树链剖分,就是赋予一个链的 ...

  4. 洛谷·【模板】点分树 | 震波【including 点分树

    初见安-这里是传送门:洛谷P6329 [模板]点分树 | 震波 一.点分树 其实你会点分治的话,点分树就是把点分治时的重心提出来重新连城一棵树. 比如当前点是u,求出子树v的重心root后将root与 ...

  5. 洛谷.4897.[模板]最小割树(Dinic)

    题目链接 最小割树模板.具体见:https://www.cnblogs.com/SovietPower/p/9734013.html. ISAP不知为啥T成0分了.. Dinic: //1566ms ...

  6. 强连通分量:洛谷P3387 模板:缩点

    传送门 顾名思义,模板awa #include <cstdio> #include <cstring> #include <cmath> #include < ...

  7. 欧拉定理(洛谷-P5091)(扩展欧拉定理实现)

    题目描述 给你三个正整数,a,m,b,你需要求:a^b mod m 输入输出格式 输入格式: 一行三个整数,a,m,b 对于全部数据: 1≤a≤10^9 1≤b≤10^{20000000} 1≤m≤1 ...

  8. 【后缀数组】洛谷P3809模板题

    题目背景 这是一道模板题. 题目描述 读入一个长度为 n n n 的由大小写英文字母或数字组成的字符串,请把这个字符串的所有非空后缀按字典序从小到大排序,然后按顺序输出后缀的第一个字符在原串中的位置. ...

  9. 洛谷 P1919 模板】A*B Problem升级版(FFT快速傅里叶)

    https://www.luogu.com.cn/problem/P1919 题目背景 本题数据已加强,请使用 FFT/NTT,不要再交 Python 代码浪费评测资源. 题目描述 给你两个正整数 a ...

  10. 洛谷 p3372 模板-线段树 1

    题目描述 如题,已知一个数列,你需要进行下面两种操作: 1.将某区间每一个数加上x 2.求出某区间每一个数的和 输入输出格式 输入格式: 第一行包含两个整数N.M,分别表示该数列数字的个数和操作的总个 ...

最新文章

  1. Openjudge2729 Blah数集(单调队列)
  2. Hibernate 基本配置文件+基本增删改查
  3. 专访 | 执拗、纯粹的网易阮良,和他的梦想团队
  4. 将计算机设置成交换机主机名,CISCO2950交换机的配置(设置密码、IP地址、主机名)...
  5. python 连 mongodb
  6. [Java基础]对象(反)序列化流
  7. win7计算机用户文件,win7系统用户文件夹改名的图文教程
  8. inventor弧度怎么标注_家里房间太大,WiFi信号覆盖不了怎么办?网件新作:分身术...
  9. 影牛社区短视频影视APP源码
  10. 行政区划编码转换区域名工具类
  11. python参数估计(一个总体均值)
  12. 手把手教你写数独计算器(1)
  13. 计算机视觉动画视频插帧难点、流程及改进
  14. 波段测试软件,超好用的波段副图(通达信公式 副图 源码 测试图)
  15. 【计算机视觉】简述对LFT-Net(大场景点云分割)的理解
  16. 100种送给老婆的生日礼物,看看有没有你需要的!
  17. 【虚拟试衣论文笔记】M3D-VTON: A Monocular-to-3D Virtual Try-On Network
  18. 独家揭秘!抖音爆款漫画变身特效的背后技术
  19. java 批量删除数据_一种批量删除数据的方法
  20. 安吉尔直饮水设备保质交付,深圳湾公园直饮水保障完成

热门文章

  1. 抽象类,接口,魔术方法
  2. Spring Boot与Docker(一):微服务架构和容器化概述
  3. 【Best Time to Buy and Sell Stock II】cpp
  4. 【转】AsyncTask的用法
  5. Windows说明Linux分区和挂载点
  6. Visual Basic 2005 中的程式語言加強功能
  7. Netty源代码学习——EventLoopGroup原理:NioEventLoopGroup分析
  8. mysql 事务、隔离级别
  9. MUI class=mui-switch开关 JQuery 控制开关
  10. eFrameWork学习笔记-eList