关于雷劈数的一些研究

upd on 2024/9/15: 拿 py 写了一下，又短又快，cpp 再见！！！

一、雷劈数的定义

背景：有个数学家走在路上看见一个 3025 的路牌被劈成 30 和 25 了，他发现 \((30+25)^2=3025\)，因此称这种数为雷劈数。

比较小的雷劈数有 \(81=(8+1)^2,100=(10+0)^2\)。

雷劈数的定义大概为：将数 \(N\) 的十进制表示从某处分成两半 \(a\) 和 \(b\)（\(b\) 可以包含前导 \(0\)，但不能为空），那么 \((a+b)^2=N\)。

二、雷劈数的求法

因为 \(b\) 可以包含前导 \(0\)，所以考虑直接枚举 \(b\) 的位数 \(n\)。这个时候十进制拼接就可以表示为 \(10^na+b\)。

于是我们的方程变成了这样：\((a+b)^2=10^na+b\)

展开移项：\(a^2+2ab+b^2-10^na-b=0\)

\(a^2+(2b-10^n)a+b^2-b=0\)

这是一个关于 \(a\) 的二次方程，想要有整数解必须满足判别式 \(\Delta=(2b-10^n)^2-4(b^2-b)\) 为完全平方数，设为 \(c^2\)。

因此 \(4b^2-4\cdot10^nb+10^{2n}-4b^2+4b=c^2\)，即 \(4(10^n-1)b=10^{2n}-c^2\)。

\(b\) 也要是整数，因此 \(10^{2n}-c^2\) 需要是 \(4(10^n-1)\) 的倍数。我们只需要找出符合这个条件的 \(c\) 即可。

首先显然 \(10^{2n}\) 是 \(4\) 的倍数，因此 \(c\) 只要是偶数就可以满足 \(4\) 的条件，我们之后再带上这个偶数的条件。现在去掉 \(4\) 之后，我们就需要找出 \(c^2\equiv10^{2n}\equiv1\pmod{10^n-1}\)，因为 \(10^n\equiv 1\pmod{10^n-1}\)。

首先 \(P=10^n-1\) 显然不是一个质数，我们先考虑将其质因数分解为 \(\prod\limits_i p_i^{a_i}\)。如果我们对于所有的 \(p_i^{a_i}\) 求出它们的合法解，那么将所有的可能组合全都用 exCRT 合并起来就可以获得 \(P\) 的所有合法解。

现在问题变成求 \(c^2\equiv 1\pmod{p_i^{a_i}}\)。注意这不是一般的二次剩余问题，因为我们寻找的数的平方是 \(1\)。

显然的，这种情况下 \(c\) 只能为 \(1\) 或 \(p_i^{a_i}-1\)，把 \(1\) 移到 \(c^2\) 处然后因式分解即可证明。

于是我们求出了所有合法的 \(c\)。如果遇到了奇数跳过即可。

对于每一个 \(c\)，其有唯一对应的两组对偶解：\(b=\dfrac{10^{2n}-c^2}{4(10^n-1)},a=\dfrac{10^n\pm c}2-b\)。

容易发现满足上面条件的 \(c\) 所对应的 \(a\) 和 \(b\) 一定合法。

于是整个问题就解决了，回顾一遍整体的流程：

1. 枚举 \(b\) 的位数 \(n\)；
2. 将 \(P=10^n-1\) 分解为 \(\prod\limits_{i=1}^k p_i^{a_i}\)；
3. 对于每一个 \(p_i^{a_i}\) 的两种选法，将它们任意组合，形成 \(2^k\) 种不同的解，每一种都用 exCRT 合并得到对应的 \(c\)；
4. 在合法的 \(c\) 中加上 \(P+1\)，然后枚举所有的 \(c\) 判断是否是偶数；
5. 对前几步得到的每一个偶数 \(c\) 都还原出对应的两组 \(a,b\)；
6. 最后将得到的所有 \((a+b)^2\) 从小到大排序。

实现有点困难，因为我带着高精度写的所以写了 17KB，目前 15 秒能算出 \(10^{50}\) 以内所有雷劈数。代码就不给了（

upd: 拿 python 重写了一下，只有不到 2kb 而且 0.5 秒能算 \(10^{60}\)，我回去加一下优化的高精度库应该还能加速，Code

实际上现在瓶颈在于质因数分解，在目前能够快速分解的范围内（对应到输入就是 65 位左右）基本都能一秒跑完了，所以可能也就这样了吧。其它的加速就又回归到了质因数分解而这并不是我想研究的部分，如果有人会更快速的分解 \(10^n-1\) 也可以在此基础上进行优化。

特别鸣谢：

1. liqingyang，刷视频的时候看到了这个问题；
2. Grammar__hbw，把倍数关系转化为了剩余问题；
3. zhouyuhang，提出了剩余问题的大体解决方案；
4. OEIS A102766，验证了程序的较小项的正确性；
5. 我自己，写+调了一晚上的史山；

另外求一个跑的比较快的高精度库，现在我用的这个太慢了。。