Description
求
\[\sum_{i-1}^{n}{i^k\cdot a^i}
\]
的值，其中 \(n\le10^{18},k\le 2\times 10^3,a\le 10^9\) 。
Solution
本题在《具体数学》一书中有详细的解法思路说明，可以在该书第二章 <和式> 中得到思路的启发。
对这个式子进行扰动。
令
\[S(k)=\sum_{i=1}^{n}{i^k}
\]
当 \(a= 1\) 时，有
\[S(k)=\sum_{i=1}^{n}{i^k} \=\sum_{i=1}^{n}(i+1)^k-(n+1)^k+1
\]
根据二项式定理
\[(x+y)^k =\sum_{i=0}^{k}{\binom{k}{i}\cdot x^iy^{k-i}}
\]
有
\[S(k)=\sum_{i=1}^{n}{\sum_{j=0}^{k}\binom{k}{j}\cdot i^j}-(n+1)^k+1\=\sum_{j=0}^{k}\binom{k}{j}{\sum_{i=1}^{n}i^j}-(n+1)^k+1\=\sum_{j=0}^{k}\binom{k}{j}{S(j)}-(n+1)^k+1
\]
由于 \(\binom{k}{k}=1\) ，则
\[\sum_{j=0}^{k-1}{\binom{k}{j}}S(j)-(n+1)^k+1=0
\]
由于\(\binom{k}{k-1}=k\) ，则
\[k\cdot S(k-1)=(n+1）^k-\sum_{j=0}^{k-2}{\binom{k}{j}S(j)-1}\S(k-1)=\frac{(n+1）^k-\sum_{j=0}^{k-2}{\binom{k}{j}S(j)-1}}{k}
\]
令 \(t=k-1\) ，有
\[S(t)=\frac{(n+1）^{t+1}-\sum_{j=0}^{t-1}{\binom{t+1}{j}S(j)-1}}{t+1}
\]
把 \(k\) 换回来，有
\[S(k)=\frac{(n+1）^{k+1}-\sum_{j=0}^{k-1}{\binom{k+1}{j}S(j)-1}}{k+1}\S(0)=n
\]
这个式子可以 $O(k^2) $求。
s同理的考虑 \(a\not= 1\) ，则令
\[S(k)=\sum_{i=1}^{n}i^ka^i
\]
可以解得
\[S(k)=\frac{(n+1)^k\cdot a^{n+1}-a\sum_{j=0}^{k-1}{\binom{k}{j}S(j)-a}}{a-1}\S(0)=\frac{a^{n+1}-a}{a-1}
\]
同理 \(O(k^2)\) 求解，该式和边界情况留给读者自证，方法类似。
Code
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define int long long
const int mod = 1e9 + 7;
const int K = 2e3 + 10;
int fc[K],fv[K],inv[K];
inline int ksm(int x,int y) {int ret=1;for(x%=mod;y;y>>=1,x=x*x%mod)if(y&1ll)ret=ret*x%mod;return ret;}
inline int C(const int x,const int y) {return x < y ? 0 : fc[x] * fv[y] % mod * fv[x - y] % mod;}
int f[K];
int n,a,k;
inline void init() {
fc[0] = 1;
for(int i = 1; i <= k + 1; i++) fc[i] = fc[i - 1] * i % mod;
fv[k + 1] = ksm(fc[k + 1],mod - 2);
for(int i = k + 1; i >= 1; i--) fv[i - 1] = fv[i] * i % mod;
inv[1] = 1;
for(int i = 2; i <= k + 1; i++) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
return ;
}
signed main() {
scanf("%lld%lld%lld",&n,&a,&k);
init();
if(a == 1) {
f[0] = n % mod;
for(int i = 1; i <= k; i++) {
int sum1 = ksm(n + 1,i + 1),sum2 = 0;
for(int j = 0; j < i; j++) (sum2 += C(i + 1,j) * f[j] % mod) %= mod;
f[i] = (sum1 - sum2 - 1 + mod) % mod * inv[i + 1] % mod;
}
} else {
int iv = ksm(a - 1,mod - 2);
f[0] = (ksm(a,n + 1) - a + mod) % mod * iv % mod;
for(int i = 1; i <= k; i++) {
int sum1 = ksm(n + 1,i) * ksm(a,n + 1) % mod,sum2 = 0;
for(int j = 0; j < i; j++) (sum2 += C(i,j) * f[j] % mod) %= mod;
sum2 = sum2 * a % mod;
f[i] = ((sum1 - sum2 + mod) % mod - a + mod) % mod * iv % mod;
}
}
cout << f[k] << endl;
return 0;
}