BZOJ2142:礼物（扩展Lucas）

最新推荐文章于 2021-09-13 17:14:16 发布

原创最新推荐文章于 2021-09-13 17:14:16 发布 · 520 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#lucas #中国剩余定理

中国剩余定理同时被 3 个专栏收录

6 篇文章

订阅专栏

Lucas定理及拓展

5 篇文章

订阅专栏

扩展欧几里得

4 篇文章

订阅专栏

本文介绍了一种求解大组合数模p(p非质数)的问题算法，通过扩展Lucas定理并结合中国剩余定理(CRT)来解决。文章详细解释了如何通过质因数分解和递归计算来求解特定模数下的组合数。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

传送门

题意：
求大组合数模p，p不是质数。

题解：
扩展lucas。
首先，将p质因数分解，得到

x \equiv a 1 (mod p k 1 1) x \equiv a 2 (mod p k 2 2) . . . x \equiv a n (mod p k n n)

$x\equiv a_1 \pmod {p_1^{k_1}}\\ x\equiv a_2 \pmod {p_2^{k_2}}\\ \\...\\x\equiv a_n \pmod {p_n^{k_n}}\\$
这个可以CRT解决。

关键是怎么求 $a_i$ 。

a i = C (n, m) = n ! m ! ( n - m ) !

$a_i=C(n,m)=\frac{n!}{m!(n-m)!}$

设 $n!=p_i^t*x,m!=p_i^{t_1}*x_1,(n-m)!=p_i^{t_2}*x_2$
可得 $C(n,m)=p_i^{t-t_1-t_2}*\frac{x}{x_1x_2}$

问题转化为求 $n!$ 有多少个 $p_i$ 以及 $n中不被p_i^t整除的积x$
首先
$t=\lfloor\frac{n}{p_i}\rfloor+\lfloor\frac{n}{p_i^2}\rfloor+\lfloor\frac{n}{p_i^3}\rfloor+...+\lfloor\frac{n}{p_i^k}\rfloor$

这个应该很好理解，每个 $p_i^k$ 只在前面 $k$ 个被计算。
可得到函数 $f(n)=f(\frac{n}{p})+\frac{n}{p}$ .
其次，考虑除去能被 $p_i$ 整除后的数列。必定每隔 $p^t$ 形成一个模意义下的循环。直接暴力求值就好了。
还有一部分是出去的数列除以 $p_i$ 后的部分。必定也构成一个阶乘，递归求解就好了。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int n;
ll Mod,m,w[10],sum,ans=1,modcnt,aa,bb;
const ll INF=0x3f3f3f3f3f3f3f3f;
struct node
{
    ll P,val;
    node(ll P=0,ll val=0):P(P),val(val){}
}mod[50];
inline ll power(ll a,ll b,ll MOD=INF)
{
    ll res=1;
    while(b)
    {
        if(b&1)res=res*a%MOD;
        a=a*a%MOD;
        b>>=1;
    }
    return res;
}
inline node calc(ll n,ll md,ll MD)
{
    ll res=1;
    if(n<md)
    {
        for(ll i=1;i<=n;i++)res=res*i%MD;
        return node(0,res);
    }
    ll cnt=n/MD;
    for(ll i=1;i<MD;i++)
    {
        if(!(i%md))continue;
        res=res*i%MD;
    }
    res=power(res,cnt,MD);
    for(ll i=cnt*MD+1;i<=n;i++)if(i%md)res=res*i%MD;
    node tmp=calc(n/md,md,MD);
    return node(tmp.P+n/md,tmp.val*res%MD);
}
inline void exgcd(ll x,ll y,ll &a,ll &b)
{
    if(!y){a=1,b=0;return;}
    exgcd(y,x%y,a,b);
    ll a2=a;
    a=b;
    b=(a2-b*(x/y));
}
inline ll inv(ll a,ll b)
{
    exgcd(a,b,aa,bb);
    aa=(aa%b+b)%b;
    return aa;
}
inline ll merge(ll &A,ll &B,ll c,ll d)
{
    exgcd(A,c,aa,bb);
    aa*=(d-B);
    aa=(aa%c+c)%c;
    B=(A*aa+B)%(A*c);
    A=A*c;
}
inline ll C(ll a,ll b)
{
    ll A=1,B=0;
    for(int i=1;i<=modcnt;i++)
    {
        ll MOD=(power(mod[i].val,mod[i].P));
        node a1=calc(a,mod[i].val,MOD);
        node b1=calc(b,mod[i].val,MOD);
        node b2=calc(a-b,mod[i].val,MOD);
        ll cnt=a1.P-b1.P-b2.P;ll val=a1.val*inv(b1.val,MOD)%MOD*inv(b2.val,MOD)%MOD;
        val=val*power(mod[i].val,cnt,MOD)%MOD;
        merge(A,B,MOD,val);
    }
    return A+B;
}
int main()
{
    scanf("%lld",&Mod);ll t=Mod;
    scanf("%lld%d",&m,&n);
    for(int i=1;i<=n;i++)scanf("%lld",&w[i]),sum+=w[i];
    if(sum>m)printf("Impossible\n");
    else
    {
        ll lim=min(1ll*100000,Mod);
        for(int i=2;i<=lim&&t!=1;i++)
        {
            if(!(t%i))
            {
                int cnt=0;
                while(!(t%i))cnt++,t/=i;
                mod[++modcnt]=node(cnt,i);
            }
        }
        for(int i=1;i<=n;i++)
        {
            (ans*=C(m,w[i]))%=Mod;
            m-=w[i];
        }
        printf("%lld\n",ans);
    }
}