新的 \(O(k+\log n)\) 做法。

考虑计算每个猴子对答案的贡献。

打个表：

1 1 2 4 8 16 32 ...

可以看出第 $ i $ 个猴子对答案的贡献是 \(i^k \times 2^{n-i-1}\)，特别地，最后一只猴子对答案的贡献是 \(n^k\)。

写成柿子：

\[n^k+\sum_{i=1}^{n-1}i^k \times 2^{n-i-1}
\]

\[n^k+2^{n-1} \times (\sum_{i=1}^{n-1}i^k \times (2^{-1})^i)
\]

我们只需要计算出 \(\sum_{i=1}^{n-1}i^k \times (2^{-1})^i\) 即可。

然后我们发现这个柿子是 CODECHEF qpolysum，然后就做完了。

还是把做法写一遍吧

qpolysum 和本题有不同之处，即 \(i\) 从 \(0\) 开始，不过没啥区别，因为你要减掉的是一个 \(0\)（

\[S(n)=\sum_{i=0}^{n-1}i^k \times m^i
\]

在本题中相当于 \(m=2^{-1}\)。

不过这个做法是猜了一个很奇怪的结论，并且做法来自校 OJ 讨论区（

我们猜 \(S(n)=m^n G(n) - G(0)\)，其中 \(G(x)\) 是一个不超过 \(k\) 次的多项式。

证明可以看这个 blog 我才不告诉你是我看不懂

然后差分一下：

\[S(n)-S(n-1)=m^n G(n) - m^{n-1} G(n-1) = (n-1)^k \times m^{n-1}
\]

\[G(n) = \frac {(n-1)^k + G(n-1)} m
\]

设 \(G(0)=x\)，那么我们在 \(n\) 为任何值的时候用 \(x\) 表示 \(G(n)\)。

因为这个多项式的次数最高为 \(k\)，而对一个 \(k\) 次多项式差分 \(k+1\) 次后为 \(0\)，所以我们把 \(G(x)\) 差分 \(k+1\) 次后得到：

\[\sum_{i=0}^{k+1} (-1)^{i+1} \binom {k+1} i G(k+1-i) = 0
\]

我们可以用 \(x\) 表示 \(G(0) \sim G(k+1)\)，然后解一个一元一次方程即可得到 \(x\)，带入可得到 \(G(1) \sim G(k+1)\) 的值。

现在我们可以使用拉格朗日插值计算 \(G(n)\) 了，答案就是 \(m^nG(n)-G(0)\)。

\(i^k\) 可以使用线性筛，所以复杂度是 \(O(k+\log n)\)。

什么？你卡我空间？\(10\) 个数组太多了？？？

实际上可以缩成 \(7\) 个 \(\rm int\) 数组和 \(1\) 个 \(\rm bool\) 数组。

我的线性筛使用的是直接纪录最小质因数而非纪录是否为质数，可以改成后者。

然后就是 \(q\) 和 \(p\) 数组在使用的时候，\(x\) 和 \(y\) 数组已经不会使用了，所以直接使用 \(x\) 和 \(y\) 代替 \(q\) 和 \(p\) 即可。

code:

#include<cstdio>
#include<cctype>
const int M=2e7+5,mod=1e9+7;
int k,n1,n2,top,x[M],y[M],idk[M],pri[M],fac[M],ifac[M];bool zhi[M];
int G[M];
inline void read(){
    char s;long long t1,t2;
    while(isdigit(s=getchar())){
        t1=n1*10ll+(s^48);n1=t1>=mod?t1%mod:t1;
        t2=n2*10ll+(s^48);n2=t2>=mod-1?t2%(mod-1):t2;
    }
}
inline int Add(const int&a,const int&b){
    return a+b>=mod?a+b-mod:a+b;
}
inline int Del(const int&a,const int&b){
    return b>a?a-b+mod:a-b;
}
inline int C(const int&n,const int&m){
    return 1ll*fac[n]*ifac[m]%mod*ifac[n-m]%mod;
}
inline int pow(int a,int b){
    int ans=1;
    for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;
    return ans;
}
inline void sieve(const int&M){
    register int i,j,x;idk[1]=1;
    for(i=2;i<=M;++i){
        if(!zhi[i])pri[++top]=i,idk[i]=pow(i,k);
        for(j=1;j<=top&&(x=i*pri[j])<=M;++j){
            idk[x]=1ll*idk[i]*idk[pri[j]]%mod;zhi[x]=true;
            if(!(i%pri[j]))break;
        }
    }
}
inline int Inter(const int&n){
    register int i,tmp,ans=0;
    x[0]=y[k+2]=1;
    for(i=1;i<=k+1;++i)x[i]=1ll*x[i-1]*Del(n,i)%mod;
    for(i=k+1;i>=1;--i)y[i]=1ll*y[i+1]*Del(n,i)%mod;
    for(i=1;i<=k+1;++i){
        if(k+1-i&1)ans=Del(ans,1ll*1ll*x[i-1]*y[i+1]%mod*G[i]%mod*ifac[i-1]%mod*ifac[k+1-i]%mod);
        else ans=Add(ans,1ll*1ll*x[i-1]*y[i+1]%mod*G[i]%mod*ifac[i-1]%mod*ifac[k+1-i]%mod);
    }
    return ans;
}
signed main(){
    register int i,X=0,Y=0;
    fac[0]=fac[1]=ifac[0]=ifac[1]=1;read();scanf("%d",&k);sieve(k+1);x[0]=1;y[0]=0;
    for(i=1;i<=k+1;++i)x[i]=Add(x[i-1],x[i-1]),y[i]=Add(y[i-1],idk[i-1]),y[i]=Add(y[i],y[i]);
    for(i=2;i<=k+1;++i)fac[i]=1ll*fac[i-1]*i%mod,ifac[i]=1ll*(mod-mod/i)*ifac[mod%i]%mod;
    for(i=2;i<=k+1;++i)ifac[i]=1ll*ifac[i-1]*ifac[i]%mod;
    for(i=0;i<=k+1;++i){
        if(i&1){
            X=Add(X,1ll*C(k+1,i)*x[k+1-i]%mod);
            Y=Add(Y,1ll*C(k+1,i)*y[k+1-i]%mod);
        }
        else{
            X=Del(X,1ll*C(k+1,i)*x[k+1-i]%mod);
            Y=Del(Y,1ll*C(k+1,i)*y[k+1-i]%mod);
        }
    }
    G[0]=mod-1ll*Y*pow(X,mod-2)%mod;
    for(i=1;i<=k+1;++i)G[i]=Add(1ll*x[i]*G[0]%mod,y[i]);
    printf("%d",Add(1ll*Del(1ll*pow(500000004,n2)*Inter(n1)%mod,G[0])*pow(2,n2-1)%mod,pow(n1,k)));
}