BZOJ 2142 礼物 神TM数论之最终版

鸣谢
http://www.myexception.cn/program/1879933.html
http://blog.csdn.net/popoqqq/article/details/39891263

题意：链接

方法：各种数论

解析：我希望我以后不要再看到数论题TAT

好，又是看了各种题解后写的一道题，大爷的题解说的好厉害但是好像并不是那么显而易见？所以不如自己来写一发，随便理解就好了。

首先呢要明确这道题要你干什么。

就是有n个东西，然后m个人，给每个人w[i]个物品，问你方案数。

随便写写就会归纳成

$c(n,sum)*c(sum,w[1])*c(sum-w[1],w[2])...$

一直乘下去就好，最终取余p。

我靠这不就是lucas定理么，这么(哔)的一道题还有啥做的。

等等样例？这个p居然不是质数！（谁告诉你是质数了！）

噢难道我可以把这道题理解为c(n,m)终极版么- -

好，步入正题，既然知道了如此坑爹的条件后我们怎么做这道题呢？

首先，我们可以把p分解质因数，变成

$p=p1^{a1}*p2^{a2}...pn^{an}$

其中，pi均为质数。

之后怎么搞呢？对的！中国剩余定理！

将每个$pi^{ai}$单独计算最终再用天朝剩余定理合并就好了。

可是组合数怎么计算呢？

c(n,m)我们把他看做三个阶乘，分别为n!,m!,(n-m)!

但是如果直接用逆元的话并不可以，因为有可能在这些阶乘中出现取模的数，这样的话就直接变为0了。所以我们要把这些数提取出来。

引用一名大牛的方法

$n!=t1*pi^{u1}$

$m!=t2*pi^{u2}$

$(n-m)!=t3*pi^{u3}$

所以t1,t2,t3我们怎么求呢？

假设19!,pi=3（popoqqq大爷举的栗子）

19!=1*2…*19

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

所以后面的123456又是个子问题可以递归求解，而这个前面的乱七八糟一堆东西我们只需要把它乘到一起就好了！并且我们提取出来了6个pi，现在继续递归下去的话，乱七八糟的东西再次乘到一起，pi继续提取，然后继续递归就搞定n！了！

现在我们只需要$t1*niyuan(t2)*niyuan(t3)*pi^{u1-u2-u3}$

这样我们就可以搞定所有的c(n,m)最终乘到一起就好了！

代码：

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define N 100010
#define mp make_pair
#define pa pair<ll,ll> 
#define fi first
#define se second
using namespace std ;
typedef long long ll;
ll p,n,m,w[N];
ll prime[N],mod[N],cnt[N],a[N];
int tot;
ll quick_my(ll a,ll b,ll M)
{
    ll ret=1;
    while(b)
    {
        if(b&1)ret=(ret*a)%M;
        a=(a*a)%M;
        b>>=1;
    }
    return ret;
}
void exgcd(ll a,ll b,ll &x,ll &y,ll &gcd)
{
    if(!b)
    {
        x=1,y=0,gcd=a;
        return ;
    }
    exgcd(b,a%b,y,x,gcd);
    y=y-a/b*x;
}
void get_factor(ll x)
{
    for(int i=2;i*i<=x;i++)
    {
        if(x%i==0)
        {
            prime[++tot]=i,cnt[tot]=0,mod[tot]=1;
            while(x%i==0){x/=i,cnt[tot]++,mod[tot]*=i;}
        }
    }
    if(x>1)prime[++tot]=x,cnt[tot]=1,mod[tot]=x;
}
ll inverse(ll a,ll b)
{
    ll xx,yy,d;
    exgcd(a,b,xx,yy,d);
    return (xx%b+b)%b; 
}
pa fact(ll k,ll n)
{
    if(n==0)return mp(0,1);
    int x=n/prime[k],y=n/mod[k];
    ll ans=1;
    if(y)
    {
        for(int i=2;i<mod[k];i++)
        {
            if(i%prime[k]!=0)ans=(ans*1ll*i)%mod[k];
        }
        ans=quick_my(ans,y,mod[k]);
    }
    for(int i=y*mod[k]+1;i<=n;i++)
    {
        if(i%prime[k]!=0)ans=ans*1ll*i%p;
    }
    pa tmp=fact(k,x);
    return mp(x+tmp.fi,ans*tmp.se%p);
}
ll calc(int k,ll n,ll m)
{
    if(n<m)return 0;
    pa a=fact(k,n),b=fact(k,m),c=fact(k,n-m);
    return quick_my(prime[k],a.fi-b.fi-c.fi,mod[k])*a.se%mod[k]*inverse(b.se,mod[k])%mod[k]*inverse(c.se,mod[k])%mod[k];
}
ll china()
{
    ll gcd,y,x=0;
    for(int i=1;i<=tot;i++)
    {
        ll r=p/mod[i];
        exgcd(mod[i],r,gcd,y,gcd);
        x=(x+r*y*a[i])%p;
    }
    return (x+p)%p;
}
ll work(ll n,ll m)
{
    for(int i=1;i<=tot;i++)a[i]=calc(i,n,m);
    return china();
}
int main()
{
    scanf("%lld%lld%lld",&p,&n,&m);
    get_factor(p);
    int sum=0;
    for(int i=1;i<=m;i++)scanf("%lld",&w[i]),sum+=w[i];
    if(sum>n){printf("Impossible\n");return 0;}
    ll ans=work(n,sum)%p;
    for(int i=1;i<=m;i++)
    {
        ans=ans*work(sum,w[i])%p;
        sum-=w[i];
    }
    printf("%lld\n",ans);
}