BZOJ 4338 BJOI2015 糖果

此题其实没什么营养就提示了 ， 简单推一个公式 ， 再上中国剩余定理即可。

$$P_{C_{m+k-1}^{m}}^{\,n}$$

其中 p <script type="math/tex" id="MathJax-Element-28">p</script>指的排列公式。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <string>
#include <vector>
#include <deque>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <algorithm>
#include <cassert>

using namespace std;

typedef long long ll;

ll  n , m , k , p , res;

void exGcd(ll a , ll b , ll& d , ll& x , ll& y)
{
    if(!b) x=1 , y=0 , d = a;
    else 
    {
        exGcd(b, a%b, d, y, x);
        y -= a/b*x;
    }
}

ll range(ll a , ll p) { return ((a%p)+p)%p; }

ll inv(ll a , ll p)
{
    ll d , x , y;
    exGcd(a, p, d, x, y);
    return range(x, p);
}

ll power(ll a , ll n)
{
    ll res = 1;
    while(n--) res*=a;
    return res;
}

struct state
{
    ll t , a , p , rp;
    state(ll p=0 , ll rp=0) { t = 0; a = 1; this->p = p; this->rp = rp; }
    state operator *(ll a)
    {
        state res = *this;
        while(a%p==0) res.t++ , a/=p;
        a%=rp;
        res.a = range(res.a*a, rp);
        return res;
    }

    state operator /(ll a)
    {
        state res = *this;
        while(a%p==0) res.t-- , a/=p;
        a%=rp;
        res.a = range(res.a*inv(a, rp), rp);
        return res;
    }
};

void solve(ll p , ll t)
{
    state now = state(p , power(p, t));
    for(ll i=1;i<=m;i++) 
    {
        now = now*(m+k-i);
        now = now/i;
    }

    ll b = range(now.a*power(now.p, now.t), now.rp) , res=1;

    for(ll i=0;i<n;i++) res = range(res*(b-i), now.rp);

    ll add = range(inv(::p/now.rp, now.rp)*(::p/now.rp), ::p);

    add = range(add*res, ::p);

    ::res = range(::res+add, ::p);
}

int main(int argc, char *argv[]) {

    cin>>n>>m>>k>>p;

    ll _p = p;
    for(ll i=2;i*i<=p;i++)if(p%i==0)
    {
        ll t = 0;
        while(_p%i==0) t++ , _p/=i;
        solve(i, t);
    }
    if(_p!=1) solve(_p, 1);
    cout<<res<<endl;

    return 0;
}