zoukankan      html  css  js  c++  java
  • [笔记] 扩展卢卡斯

    对于一般的卢卡斯定理,要求 $C_n^mspace mod space P$中的$ p $为质数;

    而扩展卢卡斯,是解决$P$不为质数时的问题,因为$P$不是质数时,很多模意义下的的除法是做不了的(没有逆元);

    首先对$P$按算术基本定理分解

    $ P = Pi p_i^{c_i} $

    对下面这个进行中国剩余定理合并,便是答案

    ${x_1} equiv C_n^m space (mod space p_1^{c_1})$

    ${x_2} equiv C_n^m space (mod space p_2^{c_2})$

    ${x_3} equiv C_n^m space (mod space p_3^{c_3})$

    $cdot cdot cdot$

    ${x_r} equiv C_n^m space (mod space p_r^{c_r})$

    然后就是求$ C_n^m space mod space p_i^{c_i} $ 的值,

    $C_n^m = frac {n!}{m!(n-m)!} $但是我们无法对他取模,因为$ p_i^{c_i} $ 与分母不一定互质

    所以我们拆一下原式:

    $ C_n^m space mod space p_i^{c_i} = frac {frac{n!}{p_i^{a}}}{frac{m!}{p_i^b}*frac{(n-m)!}{p_i^c}} space * p_i^{a-b-c} space mod space p_i^{c_i} $

    如何求a-b-c呢?

    for(R i=n;i;i/=pi) ind+=i/pi;
    for(R i=m;i;i/=pi) ind-=i/pi;
    for(R i=n-m;i;i/=pi) ind-=i/pi;

    其中$ ind=a-b-c$ (不理解手玩一下就好了qwq)

    所以现在是求$ n! 剔除掉所有pi后 space mod space p_i^{c_i}$的值

    引用网络上的例子:

    $ (22!)=(1∗2∗3∗4∗5∗6∗7∗8∗9)∗(10∗11∗12∗13∗14∗15∗16∗17∗18)∗(19∗20∗21∗22) $

    剔除3因子之后

    $ (22!)=(1∗2∗4∗5∗7∗8)∗(10∗11∗13∗14∗16∗17)∗(19∗20∗22)∗3^6∗(1∗2∗3∗4∗5∗6∗7) $

    发现按照如下分组之后前$ ⌊frac{n}{p_i^{c_i}}⌋ $组在$ mod p_i^{c_i}$后结果相同,因此只需要做一组然后快速幂即可

    对于$(19∗20∗22)$这是冗余,暴力计算即可

    对于最后的$(1∗2∗3∗4∗5∗6∗7)$递归计算即可

    3因子被剔除暂不考虑 //就是上一步的$frac{n!}{p_i^{ind}}$已经被我们扔到外面去了

    代码如下

    inline ll fac(int n,int pi,int pk) {
        if(!n) return 1; R ans=1;
        for(R i=2;i<pk;++i) if(i%pi) ans=ans*i%pk;//循环节
        ans=qpow(ans,n/pk,pk); //快速幂,即循环节的个数
        for(R i=2;i<=n%pk;++i) if(i%pi) ans=ans*i%pk;//处理最后的散块
        return ans*fac(n/pi,pi,pk)%pk; //递归求解
    }

    于是代码合起来大致长这样:

    inline ll qpow(ll a,ll b,ll p) { R ret=1; a%=p;
        for(;b;b>>=1,(a*=a)%=p) if(b&1) (ret*=a)%=p; return ret;
    }
    inline void exgcd(ll a,ll b,ll& x,ll& y) {
        if(!b) {x=1,y=0; return ;} 
        exgcd(b,a%b,y,x),y-=a/b*x;
    }
    inline ll Inv(ll n,ll p) {
        R x,y; exgcd(n,p,x,y); return (x%p+p)%p;
    }
    inline ll fac(int n,int pi,int pk) {
        if(!n) return 1; R ans=1;
        for(R i=2;i<pk;++i) if(i%pi) ans=ans*i%pk;//循环节
        ans=qpow(ans,n/pk,pk); //快速幂,即循环节的个数
        for(R i=2;i<=n%pk;++i) if(i%pi) ans=ans*i%pk;//处理最后的散块
        return ans*fac(n/pi,pi,pk)%pk; //递归求解
    }
    inline ll L(int n,int m,int pi,int pk) {
        R ind=0; for(R i=n;i;i/=pi) ind+=i/pi;
        for(R i=m;i;i/=pi) ind-=i/pi;
        for(R i=n-m;i;i/=pi) ind-=i/pi;
        R N=fac(n,pi,pk),M=fac(m,pi,pk),N_M=fac(n-m,pi,pk);//分别求每个阶乘剔除掉pi的结果
        return N*Inv(M,pk)%pk*Inv(N_M,pk)%pk*qpow(pi,ind,pk)%pk;
    }
    inline ll solve(int n,int m) { R tmp=p,ans=0;
        for(R i=2;i*i<=tmp;++i) if(tmp%i==0) {
            R pk=1; while(tmp%i==0) pk*=i,tmp/=i;
            ans=(ans+L(n,m,i,pk)*Inv(p/pk,pk)%p*p/pk%p)%p; //对P算数基本定理分解并用CRT合并答案
        } if(tmp>1) ans=(ans+L(n,m,tmp,tmp)*Inv(p/tmp,tmp)%p*p/tmp%p)%p;
        return ans;
    } 

    顺便我们可以做道例题(Luogu P2183 [国家集训队]礼物&&题解

    当然洛谷的板子也可以打一打。


    2019.05.18 

  • 相关阅读:
    下载最新Silverlight 5 Beta客户端
    oracle数据库导入导出命令!
    使用SQL Server 2008提供的表分区向导
    Microsoft Visual Studio 2010 旗舰版下载地址
    用C#创建Windows服务(Windows Services)
    Socket通信:服务端发送安全策略到flash(c#)
    Microsoft Silverlight 4 Tools for Visual Studio 2010 下载地址
    Flex打印
    .NET中三种数据类型转换的区别:(type), type.Parse, Convert类
    JQUERY 常用方法大全
  • 原文地址:https://www.cnblogs.com/Jackpei/p/10884719.html
Copyright © 2011-2022 走看看