P7976「Stoi2033」园游会题解：分形的魅力
Untitled_unrevised
·
·
题解
概述
本文介绍对于题目 P7976「Stoi2033」园游会 在现在的计算模型下，一种单次询问时间复杂度 O(\log n)，空间复杂度 O(1) 的算法。
题意简述
定义函数 \chi(n) 为：
\chi(n) = \begin{cases}
0 & n \equiv 0 \pmod 3 \\
1 & n \equiv 1 \pmod 3 \\
-1 & n \equiv 2 \pmod 3
\end{cases}
$$
G(n) = \Bigg[\sum_{k = 0}^{n}\sum_{m = k}^{n}\chi\big(\mathrm C_{m}^{k}\big)\Bigg] \bmod 1732073999
$$
其中 $\mathrm C_{m}^{k} = \dbinom{m}{k} = \dfrac{m!}{k!(m-k)!}$ 表示组合数。
也就是计算杨辉三角前 $n+1$ 行所有系数通过 $\chi$ 函数映射后的和，结果对 $1732073999$ 取模。
## 解法
任何时候先通过实例来分析总是有益无害的。首先考虑杨辉三角各系数通过 $\chi$ 函数映射后的结果（$0 \le n < 27$）：
其中 `Z` 表示 $-1$，`.` 表示 $0$。
可以发现，该系数分布具有某种分形对称性。具体地，如果已经构造出了 $3^p$ 行的三角，则将该三角形复制后放在新三角形的顶角、左上边、左下角、右上边、右下角后，再取一个系数 $\times (-1)$ 的三角放在底边的位置，其余系数填 $0$，即可得到 $3^{p+1}$ 行的三角。
因此，当 $n = 3^{p}$ 时，$G(n) = (5 - 1)G(n / 3) = 4^{p}$。记住这个结论，后面要用。
现在考虑如何用这个特性快速计算 $G(n)$。
对于 $3^p < n \le 3^{p+1}$，分两种情况讨论：
1. $3^p < n \le 2 \times 3^{p}$：代表 $n$ 的取值碰到左上边和右上边的三角形，但没有碰到左下角、右下角和底边的三角形；
顶角的三角形各系数之和为 $4^p$，左上边和右上边的三角形系数之和为 $2G(n \bmod 3^p)$，合计 $4^p + 2G(n \bmod 3^p)$，问题转化为计算 $G(n \bmod 3^p)$。
2. $2 \times 3^p < n \le 3^{p + 1}$：代表 $n$ 的取值碰到左下角、右下角和底边的三角形；
顶角的三角形各系数之和为 $4^p$，左上边和右上边的三角形系数之和为 $2 \times 4^p$，左下角、右下角和底边的三角形系数之和为 $G(n \bmod 3^p)$（别忘了底边的系数乘过 $-1$），合计 $3 \times 4^p + G(n \bmod 3^p)$，问题转化为计算 $G(n \bmod 3^p)$。
递归终点是 $n = 0$，此时 $G(n) = 1$。
将上述算法由递归转递推，从低到高枚举 $n$ 的三进制位，如果为 $1$ 则对应第 1 种情况，如果为 $2$ 则对应第 2 种情况，直接更新答案即可；如果为 $0$ 则什么也不做。
此时一步递推的时间复杂度和空间复杂度皆为 $O(1)$，枚举所有三进制位总计时间复杂度 $O(\log n)$。
参考代码：
```cpp
#include <iostream>
typedef unsigned long long u64;
const u64 P = 1732073999;
const u64 pRes[36] = {
1, 4, 16, 64, 256, 1024,
4096, 16384, 65536, 262144, 1048576, 4194304,
16777216, 67108864, 268435456, 1073741824, 830819298, 1591203193,
1168590775, 1210215102, 1376712410, 310627643, 1242510572, 1505894290,
827355163, 1577346653, 1113164615, 988510462, 489893850, 227501401,
910005604, 175874418, 703497672, 1081916689, 863518758, 1722001033
};
inline void call() {
u64 n, p = 0, res = 1;
std::cin >> n;
for(; n; n /= 3, ++p) {
switch(n % 3) {
case 0: {
break;
}
case 1: {
res = (pRes[p] + 2 * res) % P;
break;
}
case 2: {
res = (3 * pRes[p] + res) % P;
break;
}
}
}
std::cout << res << std::endl;
}
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
u64 T, maxn;
std::cin >> T >> maxn;
while(T--) call();
return 0;
}
```