• 卢卡斯定理与拓展卢卡斯


    看见神犇(yxs)的博客发现自己忘了(lucas),补一下

    卢卡斯定理

    结论:

    (p)为质数时,(C_{n}^{m}=C_{n\, mod\, p}^{m\, mod\, p}*C_{n/p}^{m/p})

    证明:

    引理:

    ((1+x)^p≡1+x^p (mod\, p))

    证明:

    根据费马小定理:((1+x)^p≡(1+x)≡1+x^p(mod\, p))

    主体证明:

    [(1+x)^n≡(1+x)^{lfloorfrac{n}{p} floor p}*(1+x)^{n\, mod\, p}(mod\, p) ]

    [(1+x)^n≡(1+x^p)^{lfloorfrac{n}{p} floor}*(1+x)^{n\, mod\, p}(mod\, p) ]

    二项式定理将两侧展开

    [sumlimits_{i=0}^{n}C_{n}^{i}x^i≡sumlimits_{i=0}^{lfloorfrac{n}{p} floor}C_{lfloorfrac{n}{p} floor}^{i}x^{pi}*sumlimits_{i=0}^{n\, mod\, p}C_{n\, mod\, p}^{j}x^j(mod\, p) ]

    (pi+j=m)时,(x)指数为(m),此时(i=lfloorfrac{m}{p} floor,j=m\, mod\, p)

    将此时(i,j)代入右侧得到

    [C_{n}^{m}=C_{lfloorfrac{n}{p} floor}^{lfloorfrac{m}{p} floor}*C_{n\, mod\, p}^{m\, mod\, p}(mod\, p) ]

    递归求解复杂度为(O(plogp)),当(n,mge p)时只能使用(lucas)

    #include<bits/stdc++.h>
    using namespace std;
    #define int long long
    inline int read()
    {
    	int x=0,f=1;
    	char ch;
    	for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
    	if(ch=='-') f=0,ch=getchar();
    	while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
    	return f?x:-x;
    }
    int T,n,m,p;
    inline int fast(int x,int k)
    {
    	int ret=1;
    	while(k)
    	{
    		if(k&1) ret=ret*x%p;
    		x=x*x%p;
    		k>>=1;
    	}                                             
    	return ret;
    }
    inline int C(int n,int m)
    {
    	if(n<m) return 0;
    	int sum1=1;
    	for(int i=n-m+1;i<=n;++i) (sum1*=i)%=p;
    	int sum2=1;
    	for(int i=1;i<=m;++i) (sum2*=i)%=p;
    	return sum1*fast(sum2,p-2)%p;
    }
    inline int lucas(int n,int m)
    {
    	if(!m) return 1;
    	return lucas(n/p,m/p)*C(n%p,m%p)%p;
    }
    signed main()
    {
    	T=read();
    	while(T--)
    	{
    		n=read(),m=read(),p=read();
    		printf("%lld
    ",lucas(n+m,m));
    	}
    return 0;
    }
    

    拓展卢卡斯

    (p)不是质数……

    数学中常见处理模数非质数套路:模数分解为质数然后(crt)

    (P)素数唯一分解(P=prodlimits_{i=1}^{k}p_i^{a_i})

    [egin{cases}C_{n}^{m}\, (mod p_1^{a_1})\C_{n}^{m}\, (mod p_2^{a_2})\……\C_{n}^{m}\, (mod p_n^{a_n})end{cases} ]

    先假设当前求的是第(k)项,即(C_{n}^{m}\, (mod\, p_k^{a_k})),后面模数用(p^k)代替

    等价于(frac{n!}{(n-m)!m!}mod\, p^k)

    不能直接求((n-m)!)(n!)的逆元,因为它们不一定和(p^k)互质,可能不存在逆元

    我们把左边含有(p)的因子提出来

    [frac{frac{n!}{p^x}}{frac{m!}{p^y}frac{(n-m)!}{p^z}}p^{x-y-z}mod\, p^k ]

    (frac{n!}{p^x})可以用(exgcd)求逆元

    先来求(frac{n!}{p^x})中的(x)

    我们知道(n)里面有(lfloorfrac{n}{p} floor)(p)的倍数,那么我们每次统计出(lfloorfrac{n}{p} floor)个贡献递归求解

    (n!)中有(g(n))(p)的因子,那么

    [g(n)=lfloorfrac{n}{p} floor+g(lfloorfrac{n}{p} floor) ]

    边界为(g(x)=0\,(x<p))

    然后考虑求(frac{n!}{p^x})的值(f(n))

    由于是在模(p^k)的意义下,所以

    [n!=(prodlimits_{i=1}^{p^k}i)^{lfloorfrac{n}{p^k} floor}prodlimits_{i=1}^{n\, mod\, p^k}i ]

    我们把含(p)的项提出来

    [n!=prodlimits_{i=1}^{lfloorfrac{n}{p} floor}(pi)(prodlimits_{i=1,i\,mod\,p e 0}^{p^k}i)^{lfloorfrac{n}{p^k} floor}prodlimits_{i=1,i\,mod\,p e 0}^{n\, mod\, p^k}i ]

    [n!=p^{lfloorfrac{n}{p} floor}(lfloorfrac{n}{p} floor)!(prodlimits_{i=1,i\,mod\,p e 0}^{p^k}i)^{lfloorfrac{n}{p^k} floor}prodlimits_{i=1,i\,mod\,p e 0}^{n\, mod\, p^k}i ]

    倒数第二个(prod)可以求一个(p^k)以内的乘积然后快速幂,最后一个(prod)直接暴力求

    至于(p^{lfloorfrac{n}{p} floor})我们应该扔掉,但是((lfloorfrac{n}{p} floor)!)内可能还有答案,递归求解

    [f(n)=f(lfloorfrac{n}{p} floor)(prodlimits_{i=1,i\,mod\,p e 0}^{p^k}i)^{lfloorfrac{n}{p^k} floor}prodlimits_{i=1,i\,mod\,p e 0}^{n\, mod\, p^k}i ]

    复杂度为(O(plogp))

    #include<bits/stdc++.h>
    using namespace std;
    #define int long long
    inline int read()
    {
    	int x=0,f=1;
    	char ch;
    	for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
    	if(ch=='-') f=0,ch=getchar();
    	while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
    	return f?x:-x;
    }
    int n,m,p;
    inline void exgcd(int &x,int &y,int a,int b)
    {
    	if(!b)
    	{
    		x=1,y=0;
    		return;
    	}
    	exgcd(y,x,b,a%b);
    	y-=a/b*x;
    }
    inline int fast(int x,int k,int p)
    {
    	int ret=1;
    	while(k)
    	{
    		if(k&1) ret=ret*x%p;
    		x=x*x%p;
    		k>>=1;
    	}
    	return ret;
    }
    inline int inv(int a,int b)
    {
    	int x,y;
    	exgcd(x,y,a,b);
    	return (x%b+b)%b;
    }
    inline int fac(int x,int a,int b)
    {
    	if(!x) return 1;
    	int ret=1;
    	for(int i=1;i<=b;++i)
    		if(i%a) (ret*=i)%=b;
    	ret=fast(ret,x/b,b);
    	for(int i=1;i<=x%b;++i)
    		if(i%a) (ret*=i)%=b;
    	return ret*fac(x/a,a,b)%b;
    }
    inline int C(int n,int m,int a,int b)
    {
    	int N=fac(n,a,b),M=fac(m,a,b),Z=fac(n-m,a,b);
    	int sum=0;
    	for(int i=n;i;i/=a) sum+=i/a;
    	for(int i=m;i;i/=a) sum-=i/a;
    	for(int i=n-m;i;i/=a) sum-=i/a;
    	return N*fast(a,sum,b)%b*inv(M,b)%b*inv(Z,b)%b;
    }
    inline int crt(int a,int m,int p)
    {
    	return inv(m/p,p)*(m/p)*a;
    }
    inline int exlucas(int n,int m,int p)
    {
    	int t=p,ret=0;
    	for(int i=2;i*i<=p;++i)
    	{
    		int k=1;
    		while(t%i==0) t/=i,k*=i;
    		if(k^1) (ret+=crt(C(n,m,i,k),p,k))%=p;
    	}
    	if(t>1) (ret+=crt(C(n,m,t,t),p,t))%=p;
    	return ret;
    }
    signed main()
    {
    	n=read(),m=read(),p=read();
    	printf("%lld
    ",exlucas(n,m,p));
    return 0;
    }
    
  • 相关阅读:
    zmodem协议的使用(SecureCRT)
    Centos6.5下安装php
    CentOS6.5下Mysql数据库的安装与配置
    YUM更换源--yum找不到安装包(转)
    百度Ueditor编辑器的Html模式自动替换样式的解决方法
    TextArea里Placeholder换行问题
    tp可用的超强第三方图表类库-JpGraph
    确保 PHP 应用程序的安全
    在线网站速度测试
    几款整站抓取的工具
  • 原文地址:https://www.cnblogs.com/knife-rose/p/12143048.html
Copyright © 2020-2023  润新知