XLkxc

Time Limit: 20 Sec Memory Limit: 128 MB

Description

给定 k,a,n,d,p
  f(i)=1^k+2^k+3^k+…+i^k
  g(x)=f(1)+f(2)+f(3)+…+f(x)
  求(g(a)+g(a+d)+g(a+2d)+…+g(a+nd))mod p

Input

第一行数据组数,(保证小于6)
  以下每行四个整数 k,a,n,d

Output

每行一个结果。

Sample Input

5
 1 1 1 1
 1 1 1 1
 1 1 1 1
 1 1 1 1
 1 1 1 1

Sample Output

5
 5
 5
 5
 5

HINT

1<=k<=123
  0<=a,n,d<=123456789
  p==1234567891

Main idea

给定k,a,n,d,求img

Solution

我们可以令img

然后推一波式子,再令img

那么显然有img

然后我们通过若干次差分,发现g在差分k+3次时全为0,那么g就是一个k+2次多项式;f在差分k+5次时全为0,那么f就是一个k+4次多项式

我们通过拉格朗日插值法插g,得到k+5个f的值,然后再插值f就可以得到答案了。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#include<bits/stdc++.h>
using namespace std;
typedef long long s64;
const int ONE=1001;
const s64 MOD=1234567891;

int T;
int k,a,n,d;
int g[ONE],f[ONE];
int inv[ONE],U[ONE],Jc[ONE];
int pre[ONE],suc[ONE];

int get()
{
int res=1,Q=1; char c;
while( (c=getchar())<48 || c>57)
if(c=='-')Q=-1;
if(Q) res=c-48;
while((c=getchar())>=48 && c<=57)
res=res*10+c-48;
return res*Q;
}

int Quickpow(int a,int b)
{
int res=1;
while(b)
{
if(b&1) res=(s64)res*a%MOD;
a=(s64)a*a%MOD;
b>>=1;
}
return res;
}

int P(int k,int i)
{
if((k-i)&1) return -1+MOD;
return 1;
}

namespace First
{
void Deal_jc(int k)
{
Jc[0]=1;
for(int i=1;i<=k;i++) Jc[i]=(s64)Jc[i-1]*i%MOD;
}

void Deal_inv(int k)
{
inv[0]=1; inv[k]=Quickpow(Jc[k],MOD-2);
for(int i=k-1;i>=1;i--) inv[i]=(s64)inv[i+1]*(i+1)%MOD;
}
}

int Final(int f[],int n,int k)
{
pre[0]=1; for(int i=1;i<=k;i++) pre[i]=(s64)pre[i-1] * (n-i+MOD) % MOD;
suc[0]=1; for(int i=1;i<=k;i++) suc[i]=(s64)suc[i-1] * (s64)(n-k+i-1+MOD) % MOD;

s64 Ans=0;
for(int i=1;i<=k;i++)
{
int Up= (s64) pre[i-1]*suc[k-i] % MOD * f[i] % MOD;
int Down= (s64) inv[i-1]*inv[k-i] % MOD;

Ans=(s64)(Ans + (s64) Up*Down % MOD * P (k,i) %MOD) % MOD;
}

return Ans;
}

int main()
{
First::Deal_jc(150); First::Deal_inv(150);
T=get();
while(T--)
{
k=get(); a=get(); n=get(); d=get();

for(int i=0;i<=k+3;i++) g[i]=Quickpow(i,k);
for(int i=1;i<=k+3;i++) g[i]=((s64)g[i-1]+g[i])%MOD;
for(int i=1;i<=k+3;i++) g[i]=((s64)g[i-1]+g[i])%MOD;
for(int i=0;i<=k+5;i++)
f[i]=((s64)f[i-1]+Final(g,(a+(s64)i*d)%MOD,k+3)) % MOD;

printf("%d\n",Final(f,n,k+5)%MOD);
}
}