P6078 [CEOI2004] Sweets

· · 题解

本题解使用 ExLucas 算法，喜提最慢解！

题意：

共 n 种糖，每种糖有 m_i 个，求选出 [l,r] 颗糖的方案数。

思路：

对于每种糖，其生成函数为（这里我们令指数表示选的个数，系数表示方案数）：

f(i) = \sum_{j=0}^{m_i} x^j

注意到：

\begin{cases} f(i) \times x = \sum_{j=1}^{m_i+1} x_j \\ f(i) - f(i) \times x = f(i) \times (1-x) =1 - x^{m_i+1}\end{cases}

得：

f(i) = \frac{1 - x^{m_i+1}}{1-x}

则最终答案式子为：

\frac{\prod\limits_{i=1}^n (1 - x^{m_i+1})}{(1-x)^n}

先考虑下面的式子：

\frac{1}{(1-x)^n} = (1-x)^{-n}

可以二项式定理：

\begin{aligned} (1-x)^{-n} &= \sum_{i=0}^{\inf} (-1)^i \times C_{-n}^i \times x^i \\ &= \sum_{i=0}^{\inf} (-1)^i \times \frac{(-n)!}{(-n-i)! \times i!} \times x^i \\ &= \sum_{i=0}^{\inf} (-1)^i \times \frac{(-n) \times (-n+1) \times \cdots \times (-n -i +1)}{i!} \times x^i \\ &=\sum_{i=0}^{\inf} (-1)^{2 \times i} \times \frac{n \times (n-1) \times \cdots \times (n+i-1)}{i!} \times x^i \\ &= \sum_{i=0}^{\inf} C_{n+i-1}^i \times x^i \end{aligned}

那么原式就转化为了：

\left( \prod\limits_{i=1}^n (1 - x^{m_i+1}) \right )\sum_{j=0}^{\inf} C_{n+j-1}^j \times x^j

那么前面的式子可以直接暴力展开计算，若得到 a \times x^b，考虑乘上后面的式子后在一个前缀 [1,m] 的贡献，即：

b \sum_{j=0}^{m-b} C_{n+j-1}^j

对于后面这个式子，根据组合数递推公式化简一下可得：

\sum_{j=0}^{m-b} C_{n+j-1}^j = C_{n+m-b}^{m-b}

对于一个区间的差分一下即可。

现在主要是对于非质数取模求出组合数，明显可以扩展卢卡斯定理。

考虑对 2004 质因数分解得 2^2 \times 3 \times 167，那么相当于求出下列式子的解：

\begin{cases} C_n^m \bmod 2^2 \\ C_n^m \bmod 3 \\ C_n^m \bmod 167\end{cases}

因为模数两两是互质的，最后可以用中国剩余定理合并，这里就不详细讲了，有需要的可以去 here 看看。

每次求组合数的常数很小，不大于为 O(p)。

时间复杂度约为 O(2^N \times p)。

完整代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
const ll M=5,mod=2004; 
inline ll read(){
    ll x=0,f=1;
    char c=getchar();
    while(c<'0'||c>'9'){
        if(c=='-')
          f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9'){
        x=(x<<1)+(x<<3)+(c^48);
        c=getchar();
    }
    return x*f;
}
inline void write(ll x){
    if(x<0){
        putchar('-');
        x=-x;
    }
    if(x>9)
      write(x/10);
    putchar(x%10+'0');
}
ll n,l,r,cnt,x,ans;
ll b[M],c[M],d[M];
map<ll,ll> H;
struct Node{
    vector<pair<ll,ll>> a;
    Node(){
        a.clear();
    }
    void init(){
        a.clear();
    }
    Node operator*(const Node&b)const{
        H.clear();
        Node ans;
        for(auto i:a)
          for(auto j:b.a)
            H[i.first+j.first]+=i.second*j.second;
        for(auto i:H)
          ans.a.emplace_back(i.first,i.second);
        sort(ans.a.begin(),ans.a.end());
        return ans;
    }
}B,A;
ll qpow(ll a,ll b,ll mod){
    ll ans=1;
    while(b){
        if(b&1ll)
          ans=(ans*a)%mod;
        a=(a*a)%mod;
        b>>=1ll;
    }
    return ans;
}
void exgcd(ll a,ll b,ll &x,ll &y){
    if(!b){
        x=1;
        y=0;
        return;
    }
    ll p;
    exgcd(b,a%b,x,y);
    p=x;
    x=y;
    y=p-(a/b)*y;
    return ;
}
ll Inv(ll a,ll mod){
    ll x=0,y=0;
    exgcd(a,mod,x,y);
    return(x%mod+mod)%mod;
}
ll solve(ll n,ll p,ll pk){
    if(!n)
      return 1;
    ll ans=1;
    for(ll i=1;i<=pk;i++)
      if(i%p)
        ans=ans*i%pk;
    ans=qpow(ans,n/pk,pk);
    for(ll i=1;i<=n%pk;i++)
      if(i%p)
        ans=ans*i%pk;
    return ans*solve(n/p,p,pk)%pk;
}
ll C(ll n,ll m,ll p,ll pk){
    if(!n||!m||n==m)
      return 1;
    if(n<m)
      return 0;
    ll X=solve(n,p,pk),Y=solve(m,p,pk),Z=solve(n-m,p,pk),cnt=0,k=n-m;
    while(n)
      n/=p,cnt+=n;
    while(m)
      m/=p,cnt-=m;
    while(k)
      k/=p,cnt-=k;
    return X*Inv(Y,pk)%pk*Inv(Z,pk)%pk*qpow(p,cnt,pk)%pk;
}
ll CRT(){
    ll M=1,ans=0;
    for(ll i=1;i<=cnt;i++)
      M*=c[i];
    for(ll i=1;i<=cnt;i++)
      d[i]=Inv(M/c[i],c[i]);
    for(ll i=1;i<=cnt;i++)
      ans+=b[i]*(M/c[i])*d[i];
    for(int i=1;i<=cnt;i++)
      b[i]=c[i]=d[i]=0;
    cnt=0;
    return (ans%M+M)%M;
}
ll ExLucas(ll n,ll m,ll p){
    if(n<m||m<0)
      return 0;
    if(!n||!m||n==m)
      return 1;
    ll t=sqrt(p);
    for(ll i=2;i<=t&&p>=1;i++){
        ll pk=1;
        while(p%i==0){
            p/=i;
            pk*=i;
        }
        if(pk>1){
            b[++cnt]=C(n,m,i,pk);
            c[cnt]=pk;
        }
    }
    if(p>1){
        b[++cnt]=C(n,m,p,p);
        c[cnt]=p;
    }
    return CRT();   
}
int main(){
    n=read(),l=read(),r=read();
    B.a.emplace_back(0,1);
    for(int i=1;i<=n;i++){
        x=read();
        A.init();
        A.a.emplace_back(0,1);
        A.a.emplace_back(x+1,-1);
        B=B*A;
    }
    for(auto v:B.a){
        x=v.first;
//      cerr<<n+r-x<<' '<<r-x<<' '<<n+l-x-1<<' '<<l-x-1<<' '<<(ExLucas(n+r-x,r-x,mod)-ExLucas(n+l-x-1,l-x-1,mod))<<'\n';
        ans=(ans+v.second*(ExLucas(n+r-x,r-x,mod)-ExLucas(n+l-x-1,l-x-1,mod))%mod+mod)%mod;
    }
    write(ans);
    return 0;
}