20240107 [Special Issue 2] 2024 元旦红包——++auto 总结

[saffahyjp]

各位久等了！鸽了一周的元旦红包活动总结终于来了！
您可以在这里回顾题目：

• 20231231 [Special Issue 1] 2024 元旦红包——++auto #3

红包领取情况

• Task 0A：32 / 100
• Task 0B：29 / 50
• Task 1A：29 / 50
• Task 1B：20 / 25
• Task 2：11 / 25
• Task 3A：9 / 100
• Task 3B：9 / 50
• Task 3C：7 / 25

Task 0

Task 0A

非常简单，唯一一处 auto 显然应当替换为某个整数，根据后面的 assert 可得答案为 20240101。

原始代码：

#include <bits/stdc++.h>
using namespace std;

int main() {
    int flag = 20240101;
    assert(flag / 10000       == 2024);
    assert(flag % 10000 / 100 ==   01);
    assert(flag         % 100 ==   01);
    cout << "auto2024OK" << flag;
    return 0;
}

Task 0B

那么答案肯定不会还是 20240101 了对吧……
注意到 C++ 中以 0 开头的整数字面量是八进制，0101 实际是十进制的 65，因此答案为 20240065。
如果你之前不知道这一点，那恭喜你从本次活动中学到了新知识！

原始代码：

#include <bits/stdc++.h>
using namespace std;

int main() {
    int x = 20240065;
    assert(x / 10000 == 2024);
    assert(x % 10000 == 0101);
    cout << "OK" << x;
    return 0;
}

Task 1

Task 1A

首先，我们容易将程序还原如下（除了一个判断条件以外）：

size_t ans = 0ull;
for (size_t i = 1ull; i <= 1000ull; i++)
    for (size_t j = 1ull; j <= 1000ull; j++)
        for (size_t k = 1ull; k <= 1000ull; k++)
            if (/* auto * auto + auto * auto == auto * auto */)
                ans++;

中间的条件应该是什么呢？我们大胆猜想这是一个判断勾股数，然后检查下答案 auto2024OK1762 的 SHA256，发现对了！
（当然，如果实在猜不出来的话，理论上可以尝试不过几千种的所有可能性……但是勾股数毕竟是最符合人类直觉的，不是吗？）

原始代码：

#include <bits/stdc++.h>
using namespace std;

int main() {
    size_t ans = 0ull;
    for (size_t i = 1ull; i <= 1000ull; i++)
        for (size_t j = 1ull; j <= 1000ull; j++)
            for (size_t k = 1ull; k <= 1000ull; k++)
                if (i * i + j * j == k * k)
                    ans++;
    cout << "auto2024OK" << ans;
    return 0;
}

Task 1B——高性能计算选手

三层循环跑 1000 还可以，但是跑 1000000 应当是不行的。不过，最内层循环显然是不必要的，可以改为开根后判断是否为 1000000 以内的整数即可。
这样就变成了 1000000 的两层循环，虽然单线程还是跑不动，但是我们可以多线程大力出奇迹！放到双路 7742 的机器上跑，很快就能出结果！

有好几位高性能计算专家都是这么直接做的。

参考代码（by @Harry-Chen）：

#include <bits/stdc++.h>
using namespace std;

int main() {
    auto count = 0ull;
    #pragma omp parallel for reduction(+:count) schedule(static)
    for (auto i = 1ull; i <= 1000000ull; i++) {
        for (auto j = 1ull; j <= 1000000ull; j++) {
            auto r = i * i + j * j;
            auto rd = (long double) r;
            auto sqd = sqrtl(rd);
            auto sqi = (uint64_t) sqd;
            if (sqi <= 1000000ull && fabsl(sqd - sqi) <= 1e-20) {
                count++;
            }
        }
    }
    cout << "auto2024OK" << count;
    return 0;
}

Task 1B——算法竞赛选手

一个经典结论是，对于本原勾股数 $i^2+j^2=k^2$（其中 $gcd(i,j,k)=1$）， $k$ 一定是奇数， $i$ 和 $j$ 中恰好有一个奇数，不妨设 $i$ 是奇数，则 $k-j$ 和 $(k-i)/2$ 都是完全平方数。该结论的证明可在中学数学范围内完成，此处略去。

其实只要知道其中两数之差是完全平方数就可以了，因此不难写出以下程序：

#include <bits/stdc++.h>
using namespace std;

void solve(size_t n) {
    size_t ans = 0;
    for (size_t k = 1ull; k <= n; k++) {
        for (size_t d = 1ull; d * d < k; d++) {
            size_t i = k - d * d;
            size_t j = sqrt((long double) (k * k - i * i)) + .5;
            if (__gcd(i, j) == 1 && i * i + j * j == k * k)
                ans += n / k * 2;
        }
    }
    cout << ans << endl;
}

int main() {
    solve(1000ull);
    solve(1000000ull);
    return 0;
}

时间复杂度 $O(n^{1.5} \log n)$，可以在几分钟内跑出 1000000。

Task 1B——算法竞赛高手

其实不用那么麻烦，直接基于 $i^2=(k+j)(k-j)$，先枚举 $i$，再枚举 $i^2$ 的两个因数，解出 $k$ 和 $j$ 并判断是否为范围内的整数即可。 由于因数个数平方的前缀和没有一个很好的界，时间复杂度不详 ，但跑 1000000 只需要几秒钟。 这个时间复杂度是因数个数平方的前缀和，是 $O(n^{1+2\log 2/\log \log n})$ 或者 $\Omega(n \log^2 n)$，by @xumingkuan，详见评论区。

#include <bits/stdc++.h>
using namespace std;

void solve(size_t n) {
    size_t ans = 0;
    vector<vector<size_t>> divisors(n + 1ull);
    for (size_t i = 1ull; i <= n; i++)
        for (size_t j = i; j <= n; j += i)
            divisors[j].push_back(i);
    vector<size_t> v(n + 1ull);
    for (size_t i = 1ull; i <= n; i++)
        for (size_t x : divisors[i])
            for (size_t y : divisors[i]) {
                if (x * y > i)
                    break;
                if (v[x * y] == i)
                    continue;
                v[x * y] = i;
                size_t diff = x * y;
                size_t sum = i * i / diff;
                size_t k = (sum + diff) / 2;
                size_t j = (sum - diff) / 2;
                if (k + j == sum && k - j == diff && 1 <= k && k <= n && 1 <= j && j <= n && i * i + j * j == k * k)
                    ++ans;
            }
    cout << ans << endl;
}

int main() {
    solve(1000ull);
    solve(1000000ull);
    return 0;
}

Task 1B——算法竞赛大佬

一位知名算法竞赛选手 @xumingkuan 给出了一个线性的做法（详见评论区）。

Task 1B——不需要写代码的做法

因为保证了这次的元旦红包不需要写代码、运行代码也能做出来，这里提供一个不需要写代码的做法：
直接去搜索“1000000 以内的”勾股数，第一条就能得到我们需要的结果。

不过我们会把 $(i,j,k)$ 和 $(j,i,k)$ 都算上，因此需要在它的答案基础上乘 2。

Task 1B——找规律选手

也可以把小规模的答案打出来（形如 0, 0, 0, 0, 2, 2, 2, 2, 2, 4, 4, 4, 6, 6, 8, 8, 10, 10, 10），在 OEIS 上搜索：

发现它提供了一个 Mathematica 的公式，直接运行一下，把 80 改成 1000000，跑一会儿就能得到答案了！
正好清华大学前段时间向学生提供了正版 Mathematica！！！

此方法由 @wangyisong1996 提供。

Task 1B——究极逃课做法

直接对着 SHA256 去爆破，对于实际上只有几百万的答案，很快就能爆破出来……

Task 2

常规做法

一看到这份代码就知道明显是上难度了，通过瞪眼法还原这个程序应该是不太可能的。
不过，注意到其中有一些故意没有被隐藏的信息，显然可以作为突破口，一个是 wmmintrin，一个是 0x01 0x02 0x04 等的 magic number。

如果去搜索这一串数，会得到满屏的 AES：

搜索 wmmintrin，也能在 Stack Overflow 上找到 AES：

具体是哪种 AES 呢？我们可以得知，这段代码是 aes-128-ecb 的加密过程，因为：

• 可以从函数中推断明文/密文是由两个 64 位整数合成的；另外 auto(&auto, &auto, 16); 也长得很像一个长度为 16 字节的 memcpy，因此宽度是 128 位；
• 顺序是 0x01 0x02 0x04……，这是加密的顺序（解密是反过来）；
• 明文/密文只有一个块，而且代码中看上去没有提供 IV 等参数，这应当是 ECB 的行为。

随便找一个网站进行 aes-128-ecb 加密即可得到答案，需要注意端序问题（明文和密钥都是 30 11 22 20 17……）。代码中明文和密钥相等，而且高 64 位和低 64 位也相等，是不是很良心？
注意，这个做法是不需要写代码或者运行代码的。

原始代码：

#include <bits/stdc++.h>
#include <wmmintrin.h>
using namespace std;

template <int rcon>
__m128i AES_128_key_exp(__m128i key) {
    __m128i keygened = _mm_aeskeygenassist_si128(key, rcon);
    keygened = _mm_shuffle_epi32(keygened, _MM_SHUFFLE(3, 3, 3, 3));
    key = _mm_xor_si128(key, _mm_slli_si128(key, 4));
    key = _mm_xor_si128(key, _mm_slli_si128(key, 4));
    key = _mm_xor_si128(key, _mm_slli_si128(key, 4));
    return _mm_xor_si128(key, keygened);
}

__uint128_t encrypt(__uint128_t key, __uint128_t in) {
    __m128i m, k;
    memcpy(&m, &in, 16);
    memcpy(&k, &key, 16);
    m = _mm_xor_si128(m, k);
    m = _mm_aesenc_si128(m, k = AES_128_key_exp<0x01>(k));
    m = _mm_aesenc_si128(m, k = AES_128_key_exp<0x02>(k));
    m = _mm_aesenc_si128(m, k = AES_128_key_exp<0x04>(k));
    m = _mm_aesenc_si128(m, k = AES_128_key_exp<0x08>(k));
    m = _mm_aesenc_si128(m, k = AES_128_key_exp<0x10>(k));
    m = _mm_aesenc_si128(m, k = AES_128_key_exp<0x20>(k));
    m = _mm_aesenc_si128(m, k = AES_128_key_exp<0x40>(k));
    m = _mm_aesenc_si128(m, k = AES_128_key_exp<0x80>(k));
    m = _mm_aesenc_si128(m, k = AES_128_key_exp<0x1b>(k));
    m = _mm_aesenclast_si128(m, k = AES_128_key_exp<0x36>(k));
    __uint128_t ret;
    memcpy(&ret, &m, 16);
    return ret;
}

__uint128_t makeBlock(size_t hi, size_t lo) {
    return ((__uint128_t) hi << 64) | lo;
}

int main() {
    cout << "autoOK" << (uint32_t) encrypt(makeBlock(0x1926081720221130ull, 0x1926081720221130ull), makeBlock(0x1926081720221130ull, 0x1926081720221130ull));
    return 0;
}

逃课做法

没想到吧，直接对着 SHA256 去爆破仍然是可行的……虽然有十亿级别的可能性，但据 @davidwangd 消息，在他的笔记本上爆破十多分钟就能出结果……

Task 3

大佬做法

一眼看出是个排序，通过两个递归函数看出是 $O(n \log^2 n)$ 的排序，通过外层递归函数看出是归并排序，满足这个特征的广为人知的排序就只剩下 Batcher's odd-even mergesort 了，输出的是它的比较次数，做完了！

常规做法

如果没有那么大佬，我们还是来看看正常做法……

首先看主函数，根据 mt19937_64 这个提示，一开始毫无疑问地是生成了一个随机序列，接下来调用了一个神秘函数，最后 assert 的是序列变成了一个有序序列，那么神秘函数就应当是一个排序的功能。

……好像分析不出什么别的有用的东西了，这个代码看上去就不是给人还原的，枚举的话工作量也太大了。
注意到 Task 3A 的存在，每个红包只值一分钱，为什么不像 Task 1 那样只留两个子任务呢？其实通过理性分析可以知道，多一个小规模的数据就能多一个验证答案的机会，Task 3A 的存在实际上是让 Task 3 整体变简单了！

那么 Task 3A 的答案怎么得到呢？我们猜想答案不会太大（考虑到 Task 3C 的 $n=2^{40}$ 的答案都不会超过 $2^{57}$，否则超过支付宝红包口令长度限制），因此直接开始爆破，得到答案是 24063。

那么，接下来的 3B 和 3C 怎么做呢？我们手里的信息只有“这是个排序”和一个小数据的答案 24063……
注意到 $n=2^{40}$ 显然不是为了能跑得出来的，肯定有办法找规律得到结果。
说到找规律，我们想到 OEIS，而且 24063 应该不是一个特别常见的数……
把 24063 和 sort 一起扔进 OEIS，我们就得到了唯一的一条答案。

$n=2^{20}$ 的答案直接就有， $n=2^{40}$ 的答案需要你用通项公式算一下，或者打开 https://oeis.org/A053545/b053545.txt 。

原始代码：

#include <bits/stdc++.h>
using namespace std;

size_t ans = 0ull;
size_t a[1024ull];

void gate(size_t x, size_t y) {
    ans++;
    if (a[x] > a[y])
        swap(a[x], a[y]);
}

void mg(size_t i, size_t x, size_t j, size_t k, size_t r) {
    if (x - i <= r && k - j <= r) {
        gate(i, j);
    } else {
        mg(i, x, j, k, r * 2ull);
        mg(i + r, x, j + r, k, r * 2ull);
        for (i += r, x -= r; i < x; i += r * 2ull)
            gate(i, i + r);
        if (i < x + r) {
            gate(i, j);
            j += r;
        }
        for (k -= r; j < k; j += r * 2ull)
            gate(j, j + r);
    }
}

void st(size_t l, size_t r) {
    if (r - l <= 1ull)
        return;
    st(l, (l + r) / 2ull);
    st((l + r) / 2ull, r);
    mg(l, (l + r) / 2ull, (l + r) / 2ull, r, 1ull);
}

int main() {
    mt19937_64 rng;
    for (size_t i = 0ull; i < 1024ull; i++)
        a[i] = rng();
    st(0ull, 1024ull);
    for (size_t i = 1ull; i < 1024ull; i++)
        assert(a[i - 1ull] <= a[i]);
    cout << "OK" << ans;
    return 0;
}

半逃课做法

先爆破出 Task 3A 和 3B，观察到答案对 $n/2$ 取模的结果都是 $n/2-1$，大力猜想 3C 也有这个性质，于是搜索空间大大减小，直接出结果！

以上所有的爆破逃课做法都是 @davidwangd 告诉我的。实在是没想到啊没想到啊！下次怎么说也把值域出大一些了（

总结

虽然也有其他的意料外的解法，但这次最没想到的还是对 SHA256 的爆破能力，甚至个人电脑爆破 32 位已经轻轻松松了……不知道这样的活动还会不会有下一次呢？如果三年以后还有的话，不知道三年以后的大家是什么样子呢……不知道三年后的电脑性能会不会让爆破变得更轻松呢

看完这篇总结，不知道你是否发现了一些之前没想到过的方法呢？或者有什么其他要发表的意见，欢迎在下面评论！

[saffahyjp]

更新了：Task 1B——找规律选手

[Harry-Chen]

我的暴力 1B 还可以砍掉一半的常数，内层 for 大于 1000000 直接 break 就好了（跑出来之后才意识到，7742 还是厉害）

[saffahyjp]

我的暴力 1B 还可以砍掉一半的常数，内层 for 大于 1000000 直接 break 就好了（跑出来之后才意识到，7742 还是厉害）

那应当也可以限制 j>i 然后 ans+=2，常数更小（
总之让我们赞美 AMD YES！

[xumingkuan]

Task 1B——算法竞赛高手

其实不用那么麻烦，直接基于 $i^2=(k+j)(k−j)$，先枚举 $i$，再枚举 $i^2$ 的两个因数，解出 $k$ 和 $j$ 并判断是否为范围内的整数即可。由于因数个数平方的前缀和没有一个很好的界，时间复杂度不详，但跑 1000000 只需要几秒钟。

因数个数是 $d(n) = O(n^{\log 2/\log \log n})$，因此因数个数平方的前缀和有一个上界是 $O(n^{1+2\log 2/\log \log n})$。
另一方面，因数个数前缀和有一个下界是 $\displaystyle\sum_{i=1}^n d(n) = \Omega(n \log n)$，因此因数个数平方的前缀和有一个下界是 $\Omega(n \log^2 n)$。
上界和下界式子直接代入 $n = 1000000$ 的话分别是 $1.47 \times 10^9$ 和 $1.91 \times 10^8$。

Task 1B——算法竞赛大佬

勾股数可以用以下公式不重不漏地生成（参考 https://en.wikipedia.org/wiki/Pythagorean_triple ）：
$$a = k \cdot (m^2 - n^2), b = k \cdot (2mn), c = k \cdot (m^2 + n^2)$$
其中 $k, m, n$ 都是正整数， $m > n$ ， $m, n$ 互质且为一奇一偶。

于是我们直接枚举 $m$ 和 $n$ 即可，勾股数组的数量是满足要求的 $k$ 的数量的两倍（因为 $a$ 和 $b$ 可互换，而且一定有 $a \neq b$ （该结论的证明可在中学数学范围内完成，此处略去））。时间复杂度 $O(n \log n)$ ，如果对gcd函数进行记忆化的话可以做到 $O(n)$ 。

参考代码：

#include <bits/stdc++.h>
using namespace std;
int main() {
	const int NUM = 1000000;
	long long ans = 0;
	for (int m = 2; m * m <= NUM; m++) {
		for (int n = m % 2 + 1; n < m; n += 2) {
			if (__gcd(m, n) == 1) {
				ans += NUM / (m * m + n * n) * 2;
			}
		}
	}
	cout << ans << endl;
	return 0;
}