2020 牛客 NOIP 赛前集训营（第四场）提高级 C – 斐波 题解

一道让我做得相当高兴的题目。

原题链接 （这是付费比赛。）

题意简述

给定长度为 $n$ 的序列 $a$。定义函数 $f(S)=\sum\limits_{T\subseteq S} F^2(\sum_{s\in T} s)$，也即“对于 $S$ 的所有子集，求 $F(\text{子集所有元素之和})\text{的平方}$ 之和”。其中 $F(x)$ 为斐波那契数列的第 $x$ 项。有$F(0)=0, F(1)=1, F(x)=F(x-1)+F(x-2) (x > 1)$。

你需要支持以下操作，操作总数为 $m$：

1 x y：将 $a_x$ 修改为 $y$；

2 l r：输出 $\sum\limits_{i=l}^{r}\sum\limits_{j=i}^{r}f(\{a_i, a_{i+1}, \cdots, a_{j}\})$，对 $998244353$ 取模。

数据范围：对于 $\text{50\%}$ 的数据，有 $n,m\leq 100$；对于 $\text{70\%}$ 的数据，有 $n,m\leq 1000$；对于 $\text{100\%}$ 的数据，有 $n,m\leq 10^5, a_i \leq 10^5$。

$\textbf{50 pts}$做法

有没有办法在多项式时间复杂度内回答单个询问呢？

枚举子集然后计算显然是 $\mathrm{O}(2^nn)$ 的，不可接受。但我们经常对类似于“序列上连续一段的子集”上的函数求和时，考虑每加入一个新的元素构成的贡献。

具体来说，假如我们求出了 $A=\sum\limits_{T \subseteq S}f(T), S=\{a_l, \cdots, a_r\}$，那么当我们考虑 $A’=\sum\limits_{T \subseteq S’}f(T), S’=\{a_l, \cdots, a_{\color{blue}r+1}\}$ 时，显然有 $$A’=A+\sum\limits_{T \subseteq S}f(T\cup \{a_{r+1}\}), S=\{a_l, \cdots, a_r\}$$。而本题中向参数集合添加一元素就等于求和，所以，若我们能够将“原来的和 $v$ ”构成的 $F^2(v)$ 的和”与“新加入的元素的贡献”独立出来，就能在 $\mathrm{O}(n)$ 的时间内计算 $\sum\limits_{i=l}^{r}f(\{a_l, \cdots, a_i\})$。

回想斐波那契数列的一些性质。现在知道非负整数 $x, \Delta$，则$$F(x+\Delta)=F(x)F(\Delta+1)+F(x-1)F(\Delta)$$。我们高兴地发现，这样一来，只需要知道 $v(T)=\sum_{s\in T} s$，就能 $\mathrm{O}(1)$ 算出加入另一个数后的 $F(\Delta+v)$ 了。那么，假如我们知道 $$\begin{aligned}& G(S)=\sum_{T \subseteq S} F(v(T))\newline & H(S)=\sum_{T \subseteq S} F(v(T)-1), S=\{a_l, \cdots, a_r\}\end{aligned}$$，记 $\Delta=a_{r+1}$，就能算出 $$\begin{aligned}G(S\cup \{\Delta\})=G(S)+G(S)F(\Delta+1)+H(S)F(\Delta) \\ H(S\cup \{\Delta\})=H(S)+G(S)F(\Delta)+H(S)F(\Delta-1)\end{aligned}$$

同理，我们扩展到平方项的情况。经过展开可以发现 $$F^2(x+\Delta)=F^2(x)F^2(\Delta+1)+2F(x)F(x-1)F(\Delta+1)F(\Delta)+F^2(x-1)F^2(\Delta)$$，其新增贡献全部可以拆分为 $F^2(v(T)),F(v(T))F(v(T-1)), F^2(v(T)-1)$ 与 $F(\Delta)$ 等项的乘积形式。具体来讲，记 $$\begin{aligned} & P(S)=\sum_{T \subseteq S} F^2(v(T))=f(S)\\ & Q(S)=\sum_{T \subseteq S} F^2(v(T)-1)\\ & K(S)=\sum_{T \subseteq S} F(v(T))F(v(T)-1), S=\{a_l,\cdots,a_r\} \end{aligned}$$，记 $a_{r+1}=\Delta$，则有 $$\begin{aligned} P(S\cup \{\Delta\})=&P(S)+P(S)F^2(\Delta+1)+Q(S)F^2(\Delta)+2K(S)F(\Delta+1)F(\Delta)\\ Q(S\cup \{\Delta\})=&Q(S)+Q(S)F^2(\Delta)+Q(S)F^2(\Delta-1)+2K(S)F(\Delta)F(\Delta-1)\\ K(S\cup \{\Delta\})=&K(S)+P(S)F(\Delta)F(\Delta+1)+Q(S)F(\Delta)F(\Delta-1)\\ &+K(S)(F^2(\Delta)+F(\Delta+1)F(\Delta-1)) \end{aligned}$$

所以，我们枚举询问区间中的 $l’\in [l, r]$，每一次遍历序列都统一计算以 $l’$ 为左端点的子区间答案之和，就可以在 $\mathrm{O}(n^2)$ 的时间内处理单个询问。

（其实评测机够快可以骗到$\textbf{70 pts}$）

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
inl ll pow2 (ll x) { return x * x % P; }
inl int calc (int l, int r) {
	// 计算f([l, l]), ..., f([l, r])的和。
	// 及时维护F^2(v)的和，F^2(v-1)的和以及F(v)F(v-1)的和，\
	分别记作F2vS, F2v_1S, FvFv_1S。
	ll F2vS = 0, F2v_1S = 1, FvFv_1S = 0, dt, res = 0,
		Fd, Fd1, Fd_1, F2vdS, F2v_1dS, FvdFv_1dS, F2d;
	for (int i = l; i <= r; ++i) {
		dt = a[i], Fd = F[dt], Fd1 = F[dt + 1],
		Fd_1 = F[dt - 1], F2d = pow2 (Fd);
		F2vdS = (F2vS * pow2 (Fd1) + F2v_1S * F2d
			 + Fd * Fd1 % P * FvFv_1S * 2) % P,
		F2v_1dS = (F2vS * F2d + F2v_1S * pow2 (Fd_1)
			 + Fd * Fd_1 % P * FvFv_1S * 2) % P,
		FvdFv_1dS = (F2vS * Fd1 % P * Fd + FvFv_1S *
			((F2d + Fd_1 * Fd1) % P) + F2v_1S * Fd_1 % P * Fd) % P;
		(F2vS += F2vdS) %= P,
		(F2v_1S += F2v_1dS) %= P,
		(FvFv_1S += FvdFv_1dS) %= P;
		res += F2vS;
	}
	return res % P;
}
inl int query (int l, int r) {
	ll res = 0;
	for (int i = l; i <= r; ++i)
		res += calc (i, r);
	return res % P;
}

$\textbf{100 pts}$做法

我们发现这个转移式子的每一项都是前一项函数值的一次项乘该项的固定系数，那就很容易使用矩阵乘法加速递推了嘛。

具体来讲，我们设状态矩阵为$\mathbf{A}=\begin{bmatrix} P(S) & Q(S) & K(S) \end{bmatrix}$，则对于$a_i=\Delta$，我们有如下转移： $$\mathbf{B}_i=\begin{bmatrix} 1+F^2(\Delta+1) & F^2(\Delta) & F(\Delta+1)F(\Delta)\\ F^2(\Delta) & 1+F^2(\Delta-1) & F(\Delta)F(\Delta-1)\\ 2F(\Delta)F(\Delta+1) & 2F(\Delta)F(\Delta-1) & 1+F^2(\Delta)+F(\Delta+1)F(\Delta-1) \end{bmatrix}$$

则 $f(\{a_l, \cdots, a_r\})=(\begin{bmatrix}0 & 1 & 0\end{bmatrix}\times \prod_{i=l}^{r}\mathbf{B}_i)_{1,1}$。同时根据矩阵乘法满足左右结合律和分配律，我们使用线段树维护区间 $[l, r]$ 的以下信息：

• 区间内所有转移矩阵的积 $\prod\limits_{i=l}^{r}\mathbf{B}_i$；
• 区间转移矩阵前缀积之和 $\sum\limits_{i=l}^{r}\prod\limits_{j=l}^{i}\mathbf{B}_j$；
• 区间转移矩阵后缀积之和（仍然是左乘）$\sum\limits_{i=l}^{r}\prod\limits_{j=i}^{r}\mathbf{B}_j$
• 区间所有子区间的转移矩阵积的和 $\sum\limits_{i=l}^{r}\sum\limits_{j=i}^{r}\prod\limits_{k=i}^{j}\mathbf{B}_k$。

在建树和查询时合并两个区间，就可利用左右子区间的信息轻松计算。单点修改同理。详细实现参见代码。

时间复杂度 $\mathrm{O}(3^3 m\log n)$，需要注意矩乘的常数问题。

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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#include <bits/stdc++.h>
using namespace std;

#define inl inline
template <typename T> inl bool read (T &x) {
	x = 0; int f = 1; char c = getchar ();
	for (; c != EOF && !isdigit (c); c = getchar ())
		if (c == '-') f = -1;
	if (c == EOF) return 0;
	for (; c != EOF && isdigit (c); c = getchar ())
		x = (x<<1) + (x<<3) + (c^48);
	x *= f; return 1;
}
template <typename T, typename... Targs>
inl bool read (T &x, Targs&... args) {
	return read (x) && read (args...); }
template <typename T> void print (T x) {
	if (x < 0) putchar ('-'), x = -x;
	if (x > 9) print (x / 10);
	putchar ('0' + x % 10);
}
template <typename T, typename... Targs>
inl void print (T x, Targs... args) {
	print (x), putchar (' '), print (args...); }

#define reint register int
#define newl putchar('\n')
typedef long long ll;
// typedef unsigned long long ull;
// typedef __int128 lll;
// typedef long double llf;
typedef pair <int, int> pint;
#define fst first
#define scd second
#define all(p) (p).begin (), (p).end ()

constexpr int N = 1e5 + 10, P = 998244353, G = N;
int n, a[N], l, r, x, type, m, F[G];

struct matrix {
	int mat[3][3], x, y;

	/* Constructors */
	using inilst = initializer_list <initializer_list<int>>;
	inl matrix () {}
	inl void init (const inilst& _mat) {
		auto it = _mat.begin ();
		for (int i = 0; i < x; ++i) {
			memset (mat[i], 0, sizeof mat[i]);
			if (it != _mat.end ()) {
				for (int j = 0; j < y; ++j)
					mat[i][j] = *(it->begin () + j);
				++it; }
		}
	}
	inl matrix (int _x, int _y, const inilst& _mat = {})
	: x (_x), y (_y) { init (_mat); }
	inl matrix (const inilst& _mat) {
		x = _mat.size (), y = 0;
		for (auto it = _mat.begin (); it != _mat.end (); ++it)
			y = max (y, (int) it->size ());
		init (_mat);
	}

	/* Operators */
	inl const int* operator [] (int _x) const { return mat[_x]; }
	inl int* operator [] (int _x) { return mat[_x]; }
	inl const matrix operator * (const matrix &b) const {
		// 特化1x3和3x3矩阵乘法以加速。
		matrix c (x, b.y);
		c.mat[0][0] = (1ll * mat[0][0] * b.mat[0][0] + 1ll * mat[0][1] * b.mat[1][0] + 1ll * mat[0][2] * b.mat[2][0]) % P,
		c.mat[0][1] = (1ll * mat[0][0] * b.mat[0][1] + 1ll * mat[0][1] * b.mat[1][1] + 1ll * mat[0][2] * b.mat[2][1]) % P,
		c.mat[0][2] = (1ll * mat[0][0] * b.mat[0][2] + 1ll * mat[0][1] * b.mat[1][2] + 1ll * mat[0][2] * b.mat[2][2]) % P;
		if (x == 3)
			c.mat[1][0] = (1ll * mat[1][0] * b.mat[0][0] + 1ll * mat[1][1] * b.mat[1][0] + 1ll * mat[1][2] * b.mat[2][0]) % P,
			c.mat[1][1] = (1ll * mat[1][0] * b.mat[0][1] + 1ll * mat[1][1] * b.mat[1][1] + 1ll * mat[1][2] * b.mat[2][1]) % P,
			c.mat[1][2] = (1ll * mat[1][0] * b.mat[0][2] + 1ll * mat[1][1] * b.mat[1][2] + 1ll * mat[1][2] * b.mat[2][2]) % P,
			c.mat[2][0] = (1ll * mat[2][0] * b.mat[0][0] + 1ll * mat[2][1] * b.mat[1][0] + 1ll * mat[2][2] * b.mat[2][0]) % P,
			c.mat[2][1] = (1ll * mat[2][0] * b.mat[0][1] + 1ll * mat[2][1] * b.mat[1][1] + 1ll * mat[2][2] * b.mat[2][1]) % P,
			c.mat[2][2] = (1ll * mat[2][0] * b.mat[0][2] + 1ll * mat[2][1] * b.mat[1][2] + 1ll * mat[2][2] * b.mat[2][2]) % P;
		return c;
	}
	inl matrix & operator *= (const matrix &b) {
		return *this = (*this) * b; }
	inl const matrix operator + (const matrix &b) const {
		return {
			{ (mat[0][0] + b.mat[0][0]) % P, (mat[0][1] + b.mat[0][1]) % P, (mat[0][2] + b.mat[0][2]) % P },
			{ (mat[1][0] + b.mat[1][0]) % P, (mat[1][1] + b.mat[1][1]) % P, (mat[1][2] + b.mat[1][2]) % P },
			{ (mat[2][0] + b.mat[2][0]) % P, (mat[2][1] + b.mat[2][1]) % P, (mat[2][2] + b.mat[2][2]) % P }
		}; }
	inl matrix & operator += (const matrix &b) {
		return *this = *this + b; }
};

inl int pow2 (ll x) { return x * x % P; }
const matrix trans (int d) {
	int Fd = F[d], Fd_1 = F[d - 1],
		Fd1 = F[d + 1], F2d = pow2 (Fd);
	return {
		{ 1 + pow2 (Fd1), F2d, 1ll * Fd1 * Fd % P },
		{ F2d, 1 + pow2 (Fd_1), 1ll * Fd * Fd_1 % P },
		{ 2ll * Fd * Fd1 % P, 2ll * Fd * Fd_1 % P,
			(1 + F2d + 1ll * Fd1 * Fd_1 % P) % P }
	};
}

struct seg_tree {
	struct info { matrix segs/* 子区间矩阵积之和 */, pres/* 前缀积之和 */, sufs/* 后缀积之和 */, prod; };
	struct node { int l, r; info dta; } t[N<<2];
	int idx[N];

	const info colle (const info &ls, const info &rs) {
		return { ls.sufs * rs.pres + ls.segs + rs.segs,
			ls.pres + ls.prod * rs.pres,
			ls.sufs * rs.prod + rs.sufs,
			ls.prod * rs.prod }; }
	inl void renew (int x, int d) {
		auto tr = trans (d);
		t[x].dta = { tr, tr, tr, tr }; }
	#define post(x) (t[x].dta = colle (t[x<<1].dta, t[x<<1|1].dta))
	void build (int x, int l, int r) {
		t[x].l = l, t[x].r = r;
		if (l == r) return idx[l] = x, renew (x, a[l]);
		int mid = t[x].l + t[x].r >> 1;
		build (x<<1, l, mid);
		build (x<<1|1, mid + 1, r);
		post (x);
	}
	inl void upd (int pos) {
		int x = idx[pos]; renew (x, a[pos]);
		while (x > 1) x >>= 1, post (x);
	}
	info query (int x, int L, int R) {
		if (t[x].l >= L && t[x].r <= R)
			return t[x].dta;
		int mid = t[x].l + t[x].r >> 1,
			sta = ((R > mid)<<1)|(L <= mid);
		if (sta < 3) return query (x<<1|(sta-1), L, R);
		return colle (query (x<<1, L, R), query (x<<1|1, L, R));
	}
} seg;

int main () {
	/* MOCK NOIP 20220618 C - 斐波 吴秋实 */
	F[1] = 1;
	for (int i = 2; i < G; ++i)
		F[i] = (F[i - 1] + F[i - 2]) % P;

	read (n, m);
	for (int i = 1; i <= n; ++i) read (a[i]);
	seg.build (1, 1, n);
	for (int i = 1; i <= m; ++i) {
		read (type);
		if (type == 1)
			read (x), read (a[x]), seg.upd (x);
		else read (l, r), print ((matrix {{ 0, 1, 0 }}
			* seg.query (1, l, r).segs)[0][0]), newl;
	}

	return 0;
}

解二

生成函数（我不会）。

解三

lty的解法。使用通项公式 $F(n)=\dfrac{1}{\sqrt{5}}\left[\left(\dfrac{1+\sqrt{5}}{2}\right)^n-\left(\dfrac{1-\sqrt{5}}{2}\right)^n\right]$ 并利用上文方法作类似合并。

但令我相当迷惑的是，$5$ 在 $p=998244353$ 下不存在二次剩余。怎么在模意义下表示 $\sqrt{5}$ 啊？！直接存储系数吗？