洛谷P4967 黑暗打击题解

题意：

有一群生物 ccj，他们在上次的星系中，发现了一群低等生物，于是想进行一波黑暗森林打击。这群低等生物即是 \(\mathsf{Hilbert}\) 鼹鼠，生活在 \(\mathsf{Hilbert}\) 星球，住在 \(\mathsf{Hilbert}\) 曲线土壤内。

这群生物决定用最傻的办法——灌水，来淹死他们。现在“高等”生物想知道，对于 \(n\) 阶的 \(\mathsf{Hilbert}\) 曲线，从上往下灌水，能淹没几个单位面积？

这是 \(1 \sim 4\) 阶的 \(\mathsf{Hilbert}\) 曲线：

\(h_1\)，如最左图所示，是一个缺上口的正方形，这个正方形的边长为 \(1\)。从\(h_2\) 开始，按照以下方法构造曲线 \(h_i\)：将 \(h_{i-1}\) 复制四份，按 \(2\times2\) 摆放。

把左上一份逆时针转 \(90^{\circ }\)，右上一份顺时针转 \(90^{\circ }\)，然后用三条单位线段将四分曲线按照左上-左下-右下-右上的顺序连接起来。如图所示，分别展示的是 \(h_2\)，\(h_3\)，\(h_4\)。加粗的线段是额外用于连接的线段。

灌水方式：

（显然这个是 \(h_3\) 的灌水面积）绿色即为无法被灌到的地方，红色为可以灌到的地方，灰色为墙，所以答案是 \(26\)，即为样例1。

一个方格有水当且仅当在它的上，左，右方格中有至少一个方格有水，最上面一层的空格都有水。

注，此题要求对 \(9223372036854775783\) 取模

输入格式：

一个整数 \(n\)，表示这个洞穴是 \(n\) 阶的 \(\mathsf{Hilbert}\) 曲线。

输出格式：

一个整数 \(ans\)，表示有 \(ans\) 个单位面积被淹没。

数据范围：

\(n \le 10^{10000}\)，数据均为手动构造，请注意常数。

这个模数是什么逆天模数啊（我还以为是 long long 自然溢出取模

总之不可否认的是，这是一个质数，~~而且是一个大的离谱的质数~~，因此记得开个**__int128**。

首先观察这个 Hillbert 曲线第 \(n\) 阶和第 \(n-1\) 阶的联系（直接搬图了）

图中黄色表示人为的含水划分线，橙色表示没有水，蓝色表示含水。

\(n=3\) 的情况如下

\(n=4\) 的情况如下

很明显，如果把 \(n=4\) 的情况从中间拆分成四个部分，那下面的两个部分就是 \(n=3\) 的情况

而上面的两个部分更像是若干个 \(n=2\) 的情况拼接在一起。于是我们考虑以下划分方式

我们把 \(a_3,b_3,c_3\) 放到 \(n=4\) 的情况里

于是我们可以得到 \(a_4,b_4,c_4\) 的递推式 \[ \begin{aligned} &a_4=2 a_3+b_3+c_3 \\[6pt]&b_4=2 a_3+2 b_3+3 c_3+1 \\[6pt]&c_4=2 c_3+1 \end{aligned} \] 进而得到 \[ \begin{aligned} &a_n=2 a_{n-1}+b_{n-1}+c_{n-1} \\[6pt]&b_n=2 a_{n-1}+2 b_{n-1}+3 c_{n-1}+1 \\[6pt]&c_n=2 c_{n-1}+1 \\[6pt]&S_n=2 a_n+2 b_n+3 c_n+1 \end{aligned} \] 其中 \(S_n\) 表示含水格子的总面积。

看上去这道题接下来就要使用矩阵快速幂了 \[ \left[\begin{array}{llll} a_2 & b_2 & c_2 & 1 \end{array}\right] \times\left[\begin{array}{llll} 2 & 2 & 0 & 0 \\ 1 & 2 & 0 & 0 \\ 1 & 3 & 2 & 0 \\ 0 & 1 & 1 & 1 \end{array}\right]^{n-2} \] 其中 \(a_2=0,~b_2=1,~c_2=1\) ，是 \(n=2\) 时的情况。

然后我们看到了逆天的数据范围：\(n \leq 10^{10000}\) ，看来需要进一步优化

首先这个 \(c_n\) 实际上没什么必要，可以直接去掉变成 \(2^{n-1} - 1\) ，得到 \[ \begin{aligned} &a_n=2 a_{n-1}+b_{n-1}+2^{n-2} - 1 \\[6pt]&b_n=2 a_{n-1}+2 b_{n-1}+2^{n-1} - 1 + 2^{n-2} + 1 \\[6pt]&S_n=2 a_n+2 b_n+2^{n} - 1 + 2^{n-1} + 1 \end{aligned} \] 同时 \(b_n\) 实际上就是 \(S_{n-1}\) ，得到 \[ \begin{aligned} &a_n = 2a_{n-1} + S_{n-2} + 2^{n-2} - 1 \\[6pt]&S_n = 2a_n + 2S_{n-1} + 2^{n} - 1 + 2^{n-1} + 1 \end{aligned} \] 上面的改一下就是 \[ a_{n+1} = 2a_n + S_{n-1} + 2^{n-1} - 1 \] 于是 \[ S_{n} - a_{n+1} = S_{n-1} + 2^n - 1 \] 变形一下 \[ S_{n-1} - S_{n-2} - 2^{n-1} + 1 = a_{n} \] 再代回去 \[ S_n = (2S_{n-1} - 2S_{n-2} - 2^n + 2) + 2S_{n-1} + 2^n - 1 + 2^{n-1} + 1 \] 整理得（这一步有点难想到） \[ S_n + 2^n = 4(S_{n-1}+2^{n-1}) -2(S_{n-2} + 2^{n-2}) \] 设 \(f_n = S_n + 2^n\) 则有 \[ f_n=4 f_{n-1}-2 f_{n-2} \] 接着实际上有两种解法

解法1：

求一下通项公式 \[ f_n=\left(\frac{3}{2}-\sqrt{2}\right)(2-\sqrt{2})^{n-1}+\left(\frac{3}{2}+\sqrt{2}\right)(2+\sqrt{2})^{n-1} \] 得到 \[ S_n=\left(\frac{3}{2}-\sqrt{2}\right)(2-\sqrt{2})^{n-1}+\left(\frac{3}{2}+\sqrt{2}\right)(2+\sqrt{2})^{n-1}-2^n \] 于是这个问题只需要求出二次剩余就结束啦！（~~本题当然存在二次剩余~~）

~~但是没写代码~~

解法2：

听说可以证明 \[ f_n \equiv f_{n \bmod (p-1)} \quad \pmod p \] 我也不知道怎么证明（~~题解写的看不懂啊~~）

但是这道题确实可以利用扩展欧拉定理优化矩阵快速幂（一般情况下是不能用的哦！）

时间复杂度 \(\mathcal{O}(T^3 \log n)\)

代码：

#include <bits/stdc++.h>
using namespace std;
#define int __int128_t
#define gc() getchar()
#define pc(a) putchar(a)
// #define INF 0x3f3f3f3f3f3f3f3f
typedef vector< vector<int> > mat;
const int mod = 9223372036854775783ll, phi = mod - 1;
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
void add(int &x, int y) { (x += y) >= mod ? x -= mod : 0; }
#define N ((int)())

mat mul(mat &A, mat &B)
{
    int n = A.size(), m = A[0].size(), p = B[0].size();
    mat C(n, vector<int>(p));
    for(int k = 0; k < m; k++)
        for(int i = 0; i < n; i++)
            for(int j = 0; j < p; j++) add(C[i][j], A[i][k] * B[k][j] % mod);
    return C;
}
mat qpow(mat A, int k)
{
    int n = max(A.size(), A[0].size());
    mat B(n, vector<int>(n));
    for(int i = 0; i < n; i++) B[i][i] = 1;
    while(k)
    {
        if(k & 1) B = mul(B, A);
        k >>= 1; A = mul(A, A); 
    }
    return B;
}
int read()
{
    char ch = gc(); int x = 0; bool ck = false;
    for(; isdigit(ch); ch = gc())
    { 
        x = x * 10 + (ch ^ 48);
        if(x >= phi) { x %= phi; ck = true;  } 
    }
    if(ck) x += phi;
    return x;
}
void print(int x)
{
    if(x < 0) { pc('-'); x = -x; }
    if(x > 9) print(x / 10);
    pc(x % 10 + '0');
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n; n = read();
    mat A(1, vector<int>(4)), B(4, vector<int>(4)); 
    A[0][0] = 0; A[0][1] = 1; A[0][2] = 1; A[0][3] = 1;
    B[0][0] = 2; B[0][1] = 2; B[0][2] = 0; B[0][3] = 0;
    B[1][0] = 1; B[1][1] = 2; B[1][2] = 0; B[1][3] = 0;
    B[2][0] = 1; B[2][1] = 3; B[2][2] = 2; B[2][3] = 0;
    B[3][0] = 0; B[3][1] = 1; B[3][2] = 1; B[3][3] = 1;
    B = qpow(B, n - 2); A = mul(A, B);
    print((A[0][0] * 2 + A[0][1] * 2 + A[0][2] * 3 + 1) % mod);
    return 0;
}