zoukankan      html  css  js  c++  java
  • 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;
    }
    
    
  • 相关阅读:
    HDU 4539郑厂长系列故事――排兵布阵(状压DP)
    HDU 2196Computer(树形DP)
    HDU 4284Travel(状压DP)
    HDU 1520Anniversary party(树型DP)
    HDU 3920Clear All of Them I(状压DP)
    HDU 3853LOOPS(简单概率DP)
    UVA 11983 Weird Advertisement(线段树求矩形并的面积)
    POJ 2886Who Gets the Most Candies?(线段树)
    POJ 2828Buy Tickets
    HDU 1394Minimum Inversion Number(线段树)
  • 原文地址:https://www.cnblogs.com/captain1/p/10382238.html
Copyright © 2011-2022 走看看