zoukankan      html  css  js  c++  java
  • 扩展 Lucas 定理

    扩展Lucas定理

    原问题

    求:

    [C^{n}_{m}quad(mod;p) ]

    其中, (n)(m) 很大,不能够直接求阶乘。(p) 大小可接受但是不保证是质数。

    学过 (Lucas) 定理后,我们知道其限制是 (p) 必须为质数。

    我们可以分解问题:

    1. 既然去了这个条件,我们可以将 (p) 进行质因数分解,模每一个 (p_i^{k_i}),然后就相当于得到一组同余方程组。

      (CRT) 解出即可。

      本人在这里因为 (CRT) 学的不太好产生了一个错误思路,导致纠结好久,总结在此

      若有人同样误入歧途,希望有些许帮助。(如果没有纠结什么奇怪想法最好先不看

    2. 那么现在问题是解:

      [C^{n}_{m}impliesfrac{m!}{(n)!(m-n)!} quad(mod;p_i^{k_i}) ]

      问题有二: 一,(m,n) 都很大,不能直接求出阶乘。

      ​ 二,(n!,(m-n)!) 内可能有 (p_i) 这个因子,无法求出逆元。

      很难解决,但是我们发现,阶乘取模一个数后呈现一个以长 (p_i^{k_i}) 循环节不断循环的模式,似乎可以只处理一个循环节,再处理余项,两者长度都不超过 (p_i^{k_i}) ,可以接受。

      另外,从其中我们可以提出所有的 (p_i) 使式子变成这样:

      [frac{frac{m!}{p_i^{a_1}}}{frac{n!}{p_i^{a_2}} imesfrac{(m-n)!}{p_i^{a_3}}} imes p_i^{a1-a2-a3} quad(mod;p_i^{k_i}) ]

      其中 (a_1-a_2-a_3) 一定是大于 (0) 的,因为 (a) 就是从 (1) 开始一定长的自然数中 (p) 这个因子的个数。数越大的区域,(p) 因子的密度越大(因为可能有 (p) 的几次方项),那么两个小数段的 (a) 值之和,一定比等长的连续数段的 (a) 值小。

      这样就避开了逆元问题。

      那么我们来具体对一个 (n!) 展开进行量化:

      [n!equiv left. egin{align} [1 imes 2 imescdots (p_i-1) imes (p_i+1) imescdots imes (p_i^{k_i}-1)] \ imes[(1+p_i^{k_i}) imes(2+p_i^{k_i}) imescdots imes(p_i^{k_i}-1+p_i^{k_i})] end{align} ight}f 循环节 m \ imescdots imes[(1+lfloorfrac{n}{p_i^{k_i}} floor) imes p_i^{k_i}) imescdots imes n]quadf余项 m\ imes(1 imes2 imesdots imeslfloorfrac{n}{p_i} floor) imes p_i^{lfloorfrac{n}{p_i} floor}quadf 阶乘项 m\ (mod;p_i^{k_i}) ]

      出现了三部分,一部分是将有 (p) 因子的项都提出一个,形成一个新的阶乘。

      第二部分是循环的,在模 (p_i^{k_i}) 意义下,所有加在后面的 (i imes p_i^{k_i}) 都可以不要。

      式子就变成这样:

      [n!equiv p_i^{lfloor frac{n}{p_i} floor}*lfloor frac{n}{p_i} floor!*(prod_{i=1;&;p_i mid i}^{p_i^{k_i}}i)^{lfloor frac{n}{p_i^{k_i}} floor}*(prod_{i=1;&;p_i mid i}^{n\%p_i^{k_i}}i)quad(mod;p_i^{k_i}) ]

      可以在函数中处理掉余项和循环节,阶乘项中, (p_i) 的幂次是要被提出去的不用管,而 (lfloor frac{n}{p_i} floor!) 项可以再递归下去,像 (Lucas) 那样,边界是输入为 (0)

      这样就完成了对阶乘的处理。

      从这里就能看出,(Exlucas) 精髓是一种 快速对模任意数意义下的阶乘进行处理 的方法。

      这样处理出三个阶乘项,并且求出下面两个的逆元,乘起来即可。

    3. 还没完,我们还没处理后面的:

      [p_i^{a1-a2-a3} ]

      可以见到,上面处理阶乘的递归函数,每一次递归可以提取出(p)({lfloorfrac{n}{p_i} floor}) 次方。

      那么我们再建立一个递归函数 (G)

      [G(x)=lfloorfrac{n}{p_i} floor+G(lfloorfrac{n}{p_i} floor) ]

      递归边界是 (n<p_i) ,这时就没有 (p_i) 可以提取了,返回 (0)

      之后快速幂求出,乘在 (2.) 的那一坨后面。

    4. 最后,把每个 (p_i) 的解用 (CRT) 组合起来。

    5. 完毕

    (Exlucas) 的重点是一种 快速对模任意数意义下的阶乘进行处理 的方法,以及 (CRT) 进行质因数分解的组合 这一思路。

    MOD

    luogu P4720

    板子题,直接贴代码了。

    (frak code)

    #include<cstdio>
    #include<algorithm>
    using namespace std;
    typedef int int_;
    #define int long long 
    
    int n,m,p;
    
    int gcd(int a,int b){
    	return b==0?a:gcd(b,a%b);
    }
    
    void exgcd(int a,int b,int &x,int &y){
    	if(b==0){
    		x=1,y=0;
    	}
    	else{
    		exgcd(b,a%b,y,x);
    		y-=(a/b)*x;
    	}
    }
    
    int get_inv(int x,int mod){
    	int d=gcd(x,mod);
    	if(d!=1) printf("err inv %lld %lld
    ",x,mod),exit(0);
    	int e,f;
    	exgcd(x,mod,e,f);
    	e=((e%mod)+mod)%mod;
    	return e;
    }
    
    int ksm(int x,int q,int mod){
    	int ret=1;
    	x%=mod;
    	while(q>0){
    		if(q&1) (ret *= x ) %= mod ;
    		(x *= x) %= mod ;
    		q>>=1;
    	}
    	return ret;
    }
    
    int fac(int x,int pi,int pik){
    	if(x == 0) return 1ll;
    	
    	int mul=1;
    	for(int i=1;i < pik;i++){
    		if(i%pi != 0) (mul *= i) %= pik;
    	}
    	
    	mul = ksm(mul,x/pik,pik);
    	
    	for(int i=(x/pik)*pik+1;i <= x;i++){
    		if(i%pi != 0) (mul *= (i%pik)) %= pik;
    	}
    	
    	(mul *= fac(x/pi,pi,pik)) %= pik;
    	
    	return mul;
    }
    
    int g(int x,int pi){
    	if(x < pi) return 0;
    	else return x/pi+g(x/pi,pi);
    }
    
    
    int exlucas(int pi,int pik){
    	return ( ( ( fac(m,pi,pik) * get_inv(fac(n,pi,pik),pik) ) %pik * get_inv(fac(m-n,pi,pik),pik) ) %pik * ksm(pi,g(m,pi)-g(n,pi)-g(m-n,pi),pik) ) %pik;
    }
    
    
    int Crt(){
    	int ans=0;
    	int mod=p;
    	int t,ai,mi,e,f;
    	for(int i=2;i*i <= mod;i++){
    		if(mod%i != 0) continue;
    		t=1;	
    		while(mod%i == 0){
    			mod/=i;
    			t*=i;
    		}
    		ai=exlucas(i,t);
    		
    		mi=p/t;
    		
    		exgcd(mi,t,e,f);
    		e=((e%t)+t)%t;
    		(e *= ai) %= p;
    		
    		(ans += (e*mi)%p) %= p ;
    		
    	}
    	if(mod != 1){
    		t=mod;
    		ai=exlucas(mod,t);
    		
    		mi=p/t;
    		
    		exgcd(mi,t,e,f);
    		e=((e%t)+t)%t;
    		(e *= ai) %= p;
    		
    		(ans += (e*mi)%p) %= p ;
    	}
    	return ans;
    }
    
    int_ main()
    {
    	//freopen("data.in","r",stdin);
    	//freopen("my.out","w",stdout); 
    	
    	scanf("%lld %lld %lld",&m,&n,&p);
    	printf("%lld",Crt());
    	return 0;
    }
    

    内容采用“知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议”进行许可。请您在转载时注明来源及链接。


    (frak by;thorn\_)

  • 相关阅读:
    学习mongodb简单安装、连接数据库、增删改查
    第三方模块glup学习
    npm 借助第三方模块nrm切换淘宝源
    nodemon 学习
    bootstrap和ie浏览器兼容性问题怎么解决?
    所得税
    债务重组
    非货币性资产交换
    政府补助
    收入 费用 和利润
  • 原文地址:https://www.cnblogs.com/thornblog/p/12235357.html
Copyright © 2011-2022 走看看