关于雷劈数的一些研究
一、雷劈数的定义
背景：有个数学家走在路上看见一个 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
n
a
+
b
。
于是我们的方程变成了这样：
(
a
+
b
)
2
=
10
n
a
+
b
展开移项：
a
2
+
2
a
b
+
b
2
−
10
n
a
−
b
=
0
a
2
+
(
2
b
−
10
n
)
a
+
b
2
−
b
=
0
这是一个关于
a
的二次方程，想要有整数解必须满足判别式
Δ
=
(
2
b
−
10
n
)
2
−
4
(
b
2
−
b
)
为完全平方数，设为
c
2
。
因此
4
b
2
−
4
⋅
10
n
b
+
10
2
n
−
4
b
2
+
4
b
=
c
2
，即
4
(
10
n
−
1
)
b
=
10
2
n
−
c
2
。
b
也要是整数，因此
10
2
n
−
c
2
需要是
4
(
10
n
−
1
)
的倍数。我们只需要找出符合这个条件的
c
即可。
首先显然
10
2
n
是
4
的倍数，因此
c
只要是偶数就可以满足
4
的条件，我们之后再带上这个偶数的条件。现在去掉
4
之后，我们就需要找出
c
2
≡
10
2
n
≡
1
(
mod
10
n
−
1
)
，因为
10
n
≡
1
(
mod
10
n
−
1
)
。
首先
P
=
10
n
−
1
显然不是一个质数，我们先考虑将其质因数分解为
∏
i
p
a
i
i
。如果我们对于所有的
p
a
i
i
求出它们的合法解，那么将所有的可能组合全都用 exCRT 合并起来就可以获得
P
的所有合法解。
现在问题变成求
c
2
≡
1
(
mod
p
a
i
i
)
。注意这不是一般的二次剩余问题，因为我们寻找的数的平方是
1
。
显然的，这种情况下
c
只能为
1
或
p
a
i
i
−
1
，把
1
移到
c
2
处然后因式分解即可证明。
于是我们求出了所有合法的
c
。如果遇到了奇数跳过即可。
对于每一个
c
，其有唯一对应的两组对偶解：
b
=
10
2
n
−
c
2
4
(
10
n
−
1
)
,
a
=
10
n
±
c
2
−
b
。
容易发现满足上面条件的
c
所对应的
a
和
b
一定合法。
于是整个问题就解决了，回顾一遍整体的流程：
枚举
b
的位数
n
；
将
P
=
10
n
−
1
分解为
k
∏
i
=
1
p
a
i
i
；
对于每一个
p
a
i
i
的两种选法，将它们任意组合，形成
2
k
种不同的解，每一种都用 exCRT 合并得到对应的
c
；
在合法的
c
中加上
P
+
1
，然后枚举所有的
c
判断是否是偶数；
对前几步得到的每一个偶数
c
都还原出对应的两组
a
,
b
；
最后将得到的所有
(
a
+
b
)
2
从小到大排序。
实现有点困难，因为我带着高精度写的所以写了 17KB，目前 15 秒能算出
10
50
以内所有雷劈数。代码就不给了（
upd: 拿 python 重写了一下，只有不到 2kb 而且 0.5 秒能算
10
60
，我回去加一下优化的高精度库应该还能加速，Code
实际上现在瓶颈在于质因数分解，在目前能够快速分解的范围内（对应到输入就是 65 位左右）基本都能一秒跑完了，所以可能也就这样了吧。其它的加速就又回归到了质因数分解而这并不是我想研究的部分，如果有人会更快速的分解
10
n
−
1
也可以在此基础上进行优化。
特别鸣谢：
liqingyang，刷视频的时候看到了这个问题；
Grammar__hbw，把倍数关系转化为了剩余问题；
zhouyuhang，提出了剩余问题的大体解决方案；
OEIS A102766，验证了程序的较小项的正确性；
我自己，写+调了一晚上的史山；