题目链接

做法 1

可以将其看成一个 $(n+1) \times (m+1)$ 的表格,其中 $f(0,m) = 0,f(n,0) = n^K$。这个递推式的组合意义就是每次可以走向量 $(0,1)$ 和 $(1,0)$,且第一次强制走 $(0,1)$ ,走到 $(n,m)$ 的方案数,这个方案数就是 $f(n,0)$ 贡献到 $f(n,m)$ 的次数。所以可以得到答案是:

$$ \sum_{i=1}^n \binom{n-i+m-1}{m-1}i^K $$

设 $m' = m-1$,首先我们会尝试考虑将 $K$ 拆斯特林数,但是 $K$ 太大了看起来无法处理斯特林数,所以我们考虑拆前面的组合数。显然前面的组合数可以写成一个 $m'$ 次多项式 $G(n)=\binom{n+m'-1}{m'}$,现在答案也就是:

$$ \sum_{i=1}^n G(n-i)i^K \tag 1 $$

考虑继续拆式子,设 $G(x) = \sum_{i=0}^{m'} g_i x^i$,带入可以得到:

$$ \sum_{i=0}^{m'} g_i \sum_{j = 1}^n (n-j)^i j^k $$

那么我们枚举 $i$,现在重点是计算后面的东西 $F(n,k_1,k_2) = \sum_{i=1}^n (n-i)^{k_1}i^{k_2}$。由于 $k_1 \leq 30$,考虑我们对其再进行二项式展开,得到:

$$ F(n,k_1,k_2) = \sum_{i=0}^{k_1}\binom{k_1}{i} n^{i} (-1)^{k_1-i} \sum_{j=1}^n j^{k_2+k_1-i} $$

后面就是个自然数幂求和了,直接拉格朗日插值即可。复杂度 $O(mK)$。但是有个严重的问题是常数太大。

一些我使用的卡常技巧:

  • 初始化拉格朗日插值的时候要求 $O(m)$ 个 $i^k$,这里可以线性筛,也提前计算好前缀和,这样可以在内存访问连续的时候计算好
  • 拉格朗日插值很多次都是重复的,本质只会有 $O(m)$ 个不同的询问,可以记忆化。
  • 拉格朗日插值要算个 $pre_i = \prod_{i}(x-i+1)$,由于询问来的是递增的,所以 $pre_i$ 不用每次都重新算
  • 最终式子中有一项贡献是 $pre_{i-1} \times inv_{i-1}$,这个可以乘到 $pre_{i-1}$ 上,减少取模常数。

注意特判 $n=0,m=0$ 就好了。

做法 2

注意到 $(1)$ 式中 $G$ 是 $m'$ 次多项式,$i^K$ 是 $K$ 次多项式,而这个形式实际上是个卷积,所以它其实是个 $K+m'$ 次多项式,再求和大概就是 $K+m'+1$ 次多项式,直接暴力求出 $f(0,m),\ldots,f(K+m'+1,m)$ 就好了。这里其实可以直接递推!比赛的时候一直以为暴力不能求。

复杂度还是 $O(Km)$,但是常数很小

my code:

#include <bits/stdc++.h>

#define fi first
#define se second
#define DB double
#define U unsigned
#define P std::pair
#define LL long long
#define LD long double
#define pb emplace_back
#define pf emplace_front
#define MP std::make_pair
#define SZ(x) ((int)x.size())
#define all(x) x.begin(),x.end()
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(int i = a;i <= b;++i)
#define ROF(i,a,b) for(int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl

const int MAXN = 2500000 + 500;
const int ha = 1e9 + 7;

inline int qpow(int a,int n=ha-2){
    int res = 1;
    while(n){
        if(n & 1) res = 1ll*res*a%ha;
        a = 1ll*a*a%ha;
        n >>= 1;
    }
    return res;
}

LL n;int m,k;
int fac[MAXN],inv[MAXN];

inline void add(int &x,int y){
    x += y-ha;x += x>>31&ha;
}

int a[60];
int pw[MAXN];
int y[MAXN];
int dp[MAXN];

inline int C(int n,int m){
    return n < 0 || m < 0 || n < m ? 0 : 1ll*fac[n]*inv[m]%ha*inv[n-m]%ha;
}

int pre[MAXN],suf[MAXN];
int pww[30][MAXN];
bool p[MAXN];
int prime[MAXN],cnt;

inline void prework(){
    FOR(i,0,29) pww[i][1] = 1;
    FOR(i,2,MAXN-1){
        if(!p[i]){
            prime[++cnt] = i;
            pww[0][i] = qpow(i,k);
            FOR(j,1,29) pww[j][i] = 1ll*pww[j-1][i]*i%ha;
        }
        for(int j = 1;j <= cnt && 1ll*i*prime[j] < MAXN;++j){
            p[i*prime[j]] = 1;
            FOR(k,0,29) pww[k][i*prime[j]] = 1ll*pww[k][i]*pww[k][prime[j]]%ha;
            if(!(i%prime[j])) break;
        }
    }
    FOR(i,2,MAXN-1) FOR(j,0,29) add(pww[j][i],pww[j][i-1]);
}

inline int gao(int m){
    if(dp[m-::k]) return dp[m-::k];
    int ans = 0;
    FOR(i,1,m+2) y[i] = pww[m-::k][i];//qpow(i,m);
//    FOR(i,2,m+2) add(y[i],y[i-1]);
    int t = n%ha;
    if(m == ::k){
        pre[0] = suf[m+3] = 1;
        FOR(i,1,m+2) pre[i] = 1ll*(t-i+ha)*pre[i-1]%ha;
        ROF(i,m+2,1) suf[i] = 1ll*(t-i+ha)*suf[i+1]%ha;
        FOR(i,1,m+2) pre[i] = 1ll*pre[i]*inv[i]%ha;
    }
    else{
        pre[m+2] = 1ll*pre[m+1]*(t-(m+2)+ha)%ha*fac[m+1]%ha*inv[m+2]%ha;
        suf[m+3] = 1;
        ROF(i,m+2,1) suf[i] = 1ll*(t-i+ha)*suf[i+1]%ha;
    }
    FOR(i,1,m+2){
        int c = 1ll*y[i]*pre[i-1]%ha*suf[i+1]%ha*inv[m+2-i]%ha;
        if((m+2-i)&1) c = ha-c;
        add(ans,c);
    }
    return dp[m-::k] = ans;
}

inline int calc(int k2){
    int ans = 0;
    FOR(i,0,k2){
        int c = 1ll*C(k2,i)*qpow(n%ha,i)%ha*gao(k+k2-i)%ha;
        if((k2-i)&1) c = ha-c;
        add(ans,c);
    }
    return ans;
}

int main(){
    fac[0] = 1;FOR(i,1,MAXN-1) fac[i] = 1ll*fac[i-1]*i%ha;
    inv[MAXN-1] = qpow(fac[MAXN-1]);ROF(i,MAXN-2,0) inv[i] = 1ll*inv[i+1]*(i+1)%ha;
    scanf("%lld%d%d",&n,&m,&k);
    prework();
    if(n == 0) return puts("0"),0;
    if(m == 0) return printf("%d\n",qpow(n%ha,k)),0;
    a[0] = 1;
    FOR(i,0,m-2){// * (x-i)
        ROF(j,59,1) a[j] = (a[j-1]-1ll*i*a[j]%ha+ha)%ha;
        a[0] = (ha-1ll*a[0]*i%ha)%ha;
    }
    FOR(i,0,m-1) a[i] = 1ll*a[i]*inv[m-1]%ha;
    int ans = 0;
//    FOR(i,1,n) add(ans,1ll*qpow(i,k)*C(n-i+m-1,m-1)%ha);
    FOR(i,0,m-1){
        int sm = 0;
        FOR(j,i,m-1) add(sm,1ll*a[j]*qpow(m-1,j-i)%ha*C(j,i)%ha);
        a[i] = sm;
    }
    FOR(i,0,m-1){
        add(ans,1ll*a[i]*calc(i)%ha);
    }
    printf("%d\n",ans);
    return 0;
}
Last modification:July 4th, 2021 at 10:55 pm
如果觉得我的文章对你有用,请随意赞赏