problem
给定 k , a , n , d k,a,n,d k,a,n,d,定义 f ( n ) = ∑ i = 1 n i k f(n)=\sum_{i=1}^ni^k f(n)=∑i=1nik, g ( n ) = ∑ i = 1 n f ( i ) g(n)=\sum_{i=1}^nf(i) g(n)=∑i=1nf(i),求:
h ( n ) = ∑ t = 0 n g ( a + t d ) h(n)=\sum_{t=0}^ng(a+td) h(n)=t=0∑ng(a+td)
答案对 1234567891 1234567891 1234567891 取模。
数据范围: 1 ≤ k ≤ 123 1\le k\le123 1≤k≤123, 0 ≤ a , n , d ≤ 123456789 0\le a,n,d\le123456789 0≤a,n,d≤123456789。
solution
- f ( n ) f(n) f(n) 是一个自然数幂和,是 ( k + 1 ) (k+1) (k+1) 次多项式。
- g ( n ) g(n) g(n) 是 f ( n ) f(n) f(n) 的前缀和,是 ( k + 2 ) (k+2) (k+2) 次多项式。
- h ( n ) h(n) h(n) 貌似有点不好算次数。
我们可以用差分,若对一个多项式差分 k k k 次得到 0 0 0,那么它是 ( k − 1 ) (k-1) (k−1) 次多项式。差分得到 h h h 是 ( k + 4 ) (k+4) (k+4) 次的。
所以直接用拉格朗日插出 g g g,再根据 g g g 插出 h h h 即可。
小 tip:其实只需要知道大致的次数即可。如果实在不能确定准确次数,可以多取几个点算,比如取 ( k + 10 ) (k+10) (k+10) 个点。
code
#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 155
#define P 1234567891
#define ll long long
using namespace std;
int k,a,n,d,up=130,f[N],g[N],h[N],inv[N],ifac[N];
int add(int x,int y) {return (ll)x+y>=P?(ll)x+y-P:x+y;}
int dec(int x,int y) {return (ll)x-y< 0?(ll)x-y+P:x-y;}
int mul(int x,int y) {return (ll)x*y%P;}
int power(int a,int b,int ans=1){
for(;b;b>>=1,a=mul(a,a))
if(b&1) ans=mul(ans,a);
return ans;
}
void prework(){
inv[1]=1,ifac[0]=ifac[1]=1;
for(int i=2;i<=up;++i) inv[i]=mul(P-P/i,inv[P-P/i*i]);
for(int i=2;i<=up;++i) ifac[i]=mul(ifac[i-1],inv[i]);
}
int lagrange(int *a,int x){
if(x<=up) return a[x];
int ans=0;
for(int i=1;i<=up;++i){
int temp=mul(a[i],power(x-i,P-2));
temp=mul(temp,mul(ifac[i-1],(up-i)&1?(P-ifac[up-i]):ifac[up-i]));
ans=add(ans,temp);
}
for(int i=1;i<=up;++i) ans=mul(ans,x-i);
return ans;
}
int main(){
int T;
scanf("%d",&T);
prework();
while(T--){
scanf("%d%d%d%d",&k,&a,&n,&d);
for(int i=1;i<=up;++i) f[i]=add(f[i-1],power(i,k));
for(int i=1;i<=up;++i) g[i]=add(g[i-1],f[i]);
h[0]=lagrange(g,a);
for(int i=1;i<=up;++i) h[i]=add(h[i-1],lagrange(g,add(a,mul(i,d))));
printf("%d\n",lagrange(h,n));
}
return 0;
}