• Lucas定理和扩展Lucas定理


    1.Lucas定理
    首先给出式子:(C_n^m\%p = C_{lfloorfrac{n}{p} floor}^{lfloorfrac{m}{p} floor} * C_{n\%p}^{m\%p}\% p),其中p为质数。
    这里给出证明……证明是我在luogu上看到的lance1ot大佬的证明,个人认为是写的很好的,在此还要做一下补充。
    首先,对于质数p,可以保证(C_p^i(1 <= i <= p-1) equiv 0(mod p)),这个比较显然,因为组合数一定是整数,而质数p的因子只有自己和1,也就是说并没有某个数能整除它,所以答案必然是p的倍数。
    根据二项式定理,((1+x)^p = C_p^0x^0 + C_p^1x^1 + C_p^2x^2+…C_p^px^p),结合上面的结论可以知道((1+x)^p equiv 1 + x^p (mod p))
    之后我们要证明lucas定理 ,我们假设(a = kp + j,b = sp + g),那我们只要证明(C_a^b = C_k^s * C_j^g\%p)即可。
    继续从二项式定理入手。((1+x)^a equiv (1+x)^{kp+j} equiv (1+x)^{kp} * (1+x)^j equiv (1+x^p)^k * (1+x)^j(mod p))
    这时候我们观察二项式展开后第b+1项,就是(C_a^bx^b),这项显然也是可以由((1+x^p)^k)中的一项和 ((1+x)^j)中的一项相乘得到的。而且是唯一的两项,就是(C_k^sx^{sp})(C_j^gx^g) ,因为(b = sp+g),那么对于((1+x^p)^k)的展开式,显然x的指数是((s+1)p)或者更高次的指数比b大,而((s-1)p)或者更低次的,后面的指数会不够用,也无法匹配成b。所以(C_a^bx^b equiv C_k^sx^{sp}*C_j^gx^g(mod p)),那么(C_a^b equiv C_k^sx^{sp}*C_j^gx^g(mod p)),即(C_n^m\%p = C_{lfloorfrac{n}{p} floor}^{lfloorfrac{m}{p} floor} * C_{n\%p}^{m\%p}\% p)

    一开始的前提条件要求了p必须是质数。否则的话(C_p^i(1 <= i <= p-1) equiv 0(mod p))无法保证,因为p很有可能被自己的因子筛掉了。比如(C_6^4\%6),答案就是3而不是0.

    2.扩展Lucas定理。
    当我们遇到p不是质数的时候怎么办呢……
    如果p能分解成几个质数的乘积,而且这些质数的指数都是1的话,可以直接套用lucas然后用CRT合并。比如SDOI2010古代猪文
    不过如果它可以分解成质数的k次幂的乘积,这样就不行了。
    于是我们有了扩展Lucas。
    首先我们把p唯一分解,假设我们现在用(p_i^k)做模数,把所有的计算出来以后还是可以用CRT合并的。
    因为(C_n^m = frac{n!}{m!(n-m)!}),所以问题变成了如何快速在模意义下算出阶乘。
    大致的方法就是,首先我们先提取出所有p的倍数,对于n,其阶乘内部有(lfloorfrac{n}{p} floor)个p的倍数,把他们全部提取出来,结果就是(p^{lfloorfrac{n}{p} floor} * lfloorfrac{n}{p} floor!),其中(lfloorfrac{n}{p} floor!)就可以递归计算。
    对于不是p的倍数的,每个(p^k)成一个循环节,在每个循环节里面是直接把乘积算出来,最后再套上(lfloorfrac{n}{p^k} floor)的指数。最后会剩余几项,那些直接暴力乘起来就行。
    举个例子。(luogu上的例子)
    假设(n = 19,p = 3,k = 2),首先我们先把p的倍数提取出来,就变成(1 * 2 * 4 * 5 * 7 * 8 * 10 * 11 * 13 * 14 * 16 * 17 * 19 * 3^6 * 6!)
    之后可以看出(1 equiv 10(mod 9)),(2 equiv 11 (mod 9)),所以原式就变成((1 * 2 * 4 * 5 * 7 * 8)^2 * 19 * 3^6 * 6!)
    最后那个单个的19就是不在循环节里面的,但是(1 equiv 19(mod 9))嘛,所以暴力计算就好了啦。

    这样递归下去计算就可以。代码实现中略微有些不同,就是对于里面每一次计算次方的(p^{lfloorfrac{n}{p} floor}),我们不在递归的时候计算,而是全部提出来,先上下消去,最后再做乘方。这个看代码实现就好。
    这里有一道板子题国家集训队 礼物,答案显然是(C_n^{w[1]}*C_{n-w[1]}^{w[2]}…),直接套扩展lucas即可。

    #include<bits/stdc++.h>
    #define rep(i,a,n) for(ll i = a;i <= n;i++)
    #define per(i,n,a) for(ll i = n;i >= a;i--)
    #define enter putchar('
    ')
    
    using namespace std;
    typedef long long ll;
    const int M = 1000005;
    const int INF = 1000000009;
    const double eps = 1e-8;
    
    ll read()
    {
       ll ans = 0,op = 1;char ch = getchar();
       while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
       while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
       return ans * op;
    }
    
    ll n,m,p,w[105],sum,ans = 1;
    ll mul(ll a,ll b,ll t)
    {
       ll res = a * b - (ll)((long double)a / t * b + eps) * t;
       return (res % t + t) % t;
    }
    ll exgcd(ll a,ll b,ll &x,ll &y)
    {
       if(!b){x = 1,y = 0;return a;}
       ll d = exgcd(b,a%b,y,x);
       y -= a / b * x;
       return d;
    }
    ll inv(ll a,ll b)
    {
       ll x,y;
       exgcd(a,b,x,y);
       return (x % b + b) % b;
    }
    ll CRT(ll b,ll t) {return mul(mul(b,inv(p/t,t),p),p/t,p);}
    ll qpow(ll a,ll b,ll t)
    {
       ll p = 1;
       while(b)
       {
          if(b & 1) p = mul(p,a,t);
          a = mul(a,a,t),b >>= 1;
       }
       return p;
    }
    
    ll fac(ll n,ll pi,ll pk)
    {
       if(!n) return 1;
       ll res = 1;
       rep(i,2,pk) if(i % pi) res *= i,res %= pk;
       res = qpow(res,n/pk,pk);
       rep(i,2,n%pk) if(i % pi) res *= i,res %= pk;
       return res * fac(n / pi,pi,pk) % pk;
    }
    
    ll C(ll n,ll m,ll pi,ll pk)
    {
       ll d = fac(n,pi,pk),d1 = fac(m,pi,pk),d2 = fac(n-m,pi,pk);
       ll k = 0;
       for(ll i = n;i;i /= pi) k += i / pi;
       for(ll i = m;i;i /= pi) k -= i / pi;
       for(ll i = n-m;i;i /= pi) k -= i / pi;
       return mul(mul(d,inv(d1,pk),pk),mul(qpow(pi,k,pk),inv(d2,pk),pk),pk);
    }
    
    ll exlucas(ll n,ll m)
    {
       ll res = 0,tmp = p,pk;
       ll lim = sqrt(p);
       rep(i,2,lim) if(!(tmp % i))
       {
          pk = 1;
          while(!(tmp%i)) pk *= i,tmp /= i;
          res += CRT(C(n,m,i,pk),pk),res %= p;
       }
       if(tmp > 1) res += CRT(C(n,m,tmp,tmp),tmp),res %= p;
       return res;
    }
    
    int main()
    {
       p = read();
       n = read(),m = read();
       rep(i,1,m) w[i] = read(),sum += w[i];
       if(sum > n) printf("Impossible
    "),exit(0);
       rep(i,1,m) ans *= exlucas(n,w[i]),ans %= p,n -= w[i];
       printf("%lld
    ",ans);
       return 0;
    }
    
    
  • 相关阅读:
    ※剑指offer系列51:二叉搜索树的第k个结点
    ※剑指offer系列50:序列化二叉树
    sqlserver添加表注释、字段注释
    3-实体数据模型与LINQ-where&OfType
    3-实体数据模型与LINQ-Select
    Jquery 在子页面上设置父页面元素的值
    开发注意事项
    函数的进阶
    文件操作的相关内容
    基本数据类型----dict
  • 原文地址:https://www.cnblogs.com/captain1/p/10382238.html
Copyright © 2020-2023  润新知