P5488 差分与前缀和

传送门:洛谷

前题提要:包含了简单的生成函数思想以及多项式乘法,是一道不可多得的多项式好题.故记录一下.

题意:给定一个长为 n 的序列 a，求出其 k 阶差分或前缀和。结果的每一项都需要对 1004535809取模。

对于差分和前缀和我们分开来讨论.

先讨论前缀和部分:

先列出$A$序列:$(a_1,a_2,...,a_n)$,再列出前缀和$S$序列:$(s_1,s_2,...,s_n)$,为了等会的方便起见,我们将这两个序列扩展一下,分别加上一个$a_0,s_0$.(并且这两个值规定为0)
存在$s_0=0,s_1=a_0+a_1,s_n=a_0+a_1+...+a_n$
我们将其往多项式那边靠,不难发现,如果存在另一个序列$B$为$(1,1,\mathrm{1...1})$,那么此时我们的$S=A\ast B$,也就是$S$序列是$A$与$B$的卷积系数结果.
所以一阶前缀和$S1=A1\ast B1$,那么二阶$S2=S1\ast B1$,因为多项式卷积满足结合律$k$阶即为$Sk=S1\ast B{1}^k$.
此时如果多项式的科技比较强,直接套一个多项式快速幂就结束了.
但是显然我的科技树比较差,所以来讨论一下低科技的解法
因为我们无法进行多项式快速幂,所以我们考虑使用生成函数的方法将多项式转为一个函数然后在做考虑.那么此时我们的$B$就是$1+x+x^2+...+x^n$,求一下和就是$(1-x^n)/(1-x)$,为了方便起见,我们此时不妨选择$x\in (0,1)$,那么此时我们的$B$就是$1/(1-x)$,那么此时$B^k=(1-x{)}^{-k}$.
考虑使用扩展二项式定理:$$(a+b{)}^n=a^n\ast (1+\frac{b}{a}{)}^n=a^n\ast \sum _{i=0}^{\infty }\left(\genfrac{}{}{0ex}{}{n}{i}\right)\ast (\frac{b}{a}{)}^i$$
那么对于上述式子来说就是$$B^k=\sum _{i=0}^{\infty }\left(\genfrac{}{}{0ex}{}{-k}{i}\right)\ast (-1{)}^i\ast x^i$$$$\left(\genfrac{}{}{0ex}{}{-k}{i}\right)=\frac{(-k)\ast (-k-1)\ast (-k-2)\ast ...\ast (-k-i+1)}{i!}$$$$=(-1{)}^i\ast \frac{k\ast (k+1)\ast (k+2)\ast ...(k+i-1)}{i!}$$$$=(-1{)}^i\ast C_{k+i-1}^i$$

对于$C_{k+i-1}^i$,显然我们是递推来求和的(注意本题的$k$很大,只能递推来求和,无法预处理),$C_{k+i-1}^i=C_{k+i-2}^{i-1}\ast (\frac{k+i-1}{i})$.

那么此时我们的k阶前缀和就到此结束了,只要使用$NTT$解决一下多项式乘法即可.

---

下面讲一下这道题的k阶差分部分.其实大体上的解决方案和前缀和是一样的.
同样考虑得出B序列为$(1,-1,0,0,\mathrm{0...0})$ [推导方式和上述一样,此处就不在赘述了].
使用生成函数将$B^k$序列转化一下就是$(1-x{)}^k$,使用二项式定理展开就是${\sum }_{i=0}^{\infty }{C}_k^i\ast (-x{)}^k$,同样递推一下即可.

最后也是拿$NTT$做一下多项式乘法即可解决.

---

下面是具体的代码部分:

#include 
using namespace std;
typedef long long ll;
#define root 1,n,1
#define ls rt<<1
#define rs rt<<1|1
#define lson l,mid,rt<<1
#define rson mid+1,r,rt<<1|1
inline ll read() {
	ll x=0,w=1;char ch=getchar();
	for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') w=-1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
	return x*w;
}
inline void print(__int128 x){
	if(x<0) {putchar('-');x=-x;}
	if(x>9) print(x/10);
	putchar(x%10+'0');
}
#define maxn 1000000
#define int long long
const int mod=1004535809;
const double eps=1e-8;
#define	int_INF 0x3f3f3f3f
#define ll_INF 0x3f3f3f3f3f3f3f3f
int f[maxn],g1[maxn],g2[maxn];
inline ll get_read() {
	ll x=0,w=1;char ch=getchar();
	for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') w=-1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x*10%mod+ch-'0')%mod;
	return x*w;
}
int qpow(int a,int b) {
	int ans=1;
	while(b) {
		if(b&1) ans=ans*a%mod;
		b>>=1;
		a=a*a%mod;
	}
	return ans;
}
int rev[maxn];
void NTT(int *a,int n,int inv) {
	for(int i=0;i<=n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int len=1;len<=(n>>1);len<<=1) {
		int gn=qpow(inv==1?3:qpow(3,mod-2),(mod-1)/(len<<1));
		for(int i=0;i<=n;i+=(len<<1)) {
			int g0=1;
			for(int j=0;j<=len-1;j++) {
				int x=a[i+j],y=a[i+j+len]*g0%mod;
				a[i+j]=(x+y)%mod;a[i+j+len]=((x-y)%mod+mod)%mod;
				g0=g0*gn%mod;
			}
		}
	}
}
signed main() {
	int n=read();int k=get_read();int t=read();
	for(int i=1;i<=n;i++) {
		f[i]=read();
	}
	int limit=1,len=0;
	while(limit<=2*n) limit<<=1,len++;
	for(int i=0;i<=limit;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	g1[0]=1;
	for(int i=1;i<=n;i++) {
		g1[i]=g1[i-1]*(k+i-1)%mod*qpow(i,mod-2)%mod;
	}
	g2[0]=1;
	for(int i=1;i<=n;i++) {
		g2[i]=(g2[i-1]*(k-i+1)%mod*qpow(i,mod-2)%mod*-1+mod)%mod;
	}
	if(t==0) {
		NTT(f,limit,1);NTT(g1,limit,1);
		for(int i=0;i<=limit;i++) {
			f[i]=f[i]*g1[i]%mod;
		}
		NTT(f,limit,-1);
		int inv=qpow(limit,mod-2);
		for(int i=0;i<=limit;i++) {
			f[i]=f[i]*inv%mod;
		}
	}
	else {
		NTT(f,limit,1);NTT(g2,limit,1);
		for(int i=0;i<=limit;i++) {
			f[i]=f[i]*g2[i]%mod;
		}
		NTT(f,limit,-1);
		int inv=qpow(limit,mod-2);
		for(int i=0;i<=limit;i++) {
			f[i]=f[i]*inv%mod;
		}
	}
	for(int i=1;i<=n;i++) {
		cout<<f[i]<<" ";
	}
	return 0;	
}
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
91
92
93
94
95
96
97
98
99