雷劈数及其相关研究

luyuanhong

雷劈数及其相关研究

原创 陈雨翔 SDSZ CS社 2024 年 10 月 12 日 22:58 北京

一、雷劈数的定义

背景：有个数学家走在路上看见一个 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 。这个时候十进制拼接就可以表示为 a×10^n+b 。

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

展开移项：a^2+2ab+b^2-a×10^n-b=0 ，

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

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

因此 4b^2-4b×10^n+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≡10^(2n)≡1 (mod 10^n-1) ，因为 10^n≡1 (mod 10^n-1) 。

首先 p=10^n-1 显然不是一个质数，我们先考虑将其质因数分解为 ∏(i=1,m)pi^ai 。如果我们对于所有的 pi^ai 求出它们的合法解，那么将所有的可能组合全都用 exCRT 合并起来就可以获得 p 的所有合法解。

现在问题变成求 c^2≡1 (mod pi^ai) 。注意这不是一般的二次剩余问题，因为我们寻找的数的平方是 1 。

显然的，这种情况下 c 只能为 1 或 pi^ai-1 ，把 1 移到 c^2 处然后因式分解即可证明。

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

对于每一个 c ，其有唯一对应的两组对偶解：b=[10^(2n)-c^2]／[4(10^n-1)] ， a=(10^n±c)／2-b 。   

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

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

1. 枚举 b 的位数 n ；

2. 将 p=10^n-1 分解为 ∏(i=1,m)pi^ai ；

3. 对于每一个 pi^ai 的两种选法，将它们任意组合，形成 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 ，我回去加一下优化的高精度库应该还能加速。

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

SDSZ CS社

王守恩

雷劈数——卡普利加数——Kaprekar 数。

{1, 9, 45, 55, 99, 297, 703, 999, 2223, 2728, 4879, 4950, 5050, 5292, 7272, 7777, 9999, 17344, 22222, 38962, 77778, 82656, 95121, 99999, 142857, 148149, 181819, 187110, 208495, 318682, 329967,
351352, 356643, 390313, 461539, 466830, 499500, 500500, 533170, 538461, 609687, 627615, 643357, 648648, 670033, 681318, 791505, 812890, 818181, 851851, 857143, 961038, 994708, 999999,
4444444, 4927941, 5072059, 5479453,5555556, 8161912, 9372385, 9999999,11111112, 13641364,16590564, 19273023, 19773073, 24752475, 25252525, 30884184, 36363636, 38883889, 44363341,
44525548, 49995000, 50005000, 55474452, 55636659, 61116111, 63636364, 69115816, 74747475, 75247525, 80226927, 80726977, 83409436, 86358636, 88888888, 91838088, 94520547,99999999,

1. k[a_] := Module[{n = a^2}, MemberQ[Plus @@ # & /@ Select[Table[{Floor[n/10^j], FractionalPart[n/10^j] 10^j}, {j, IntegerLength[n] - 1}], #[[2]] != 0 &], a]]; Select[Range[10^9], k]

复制代码

详见——OEIS——A053816。

[code]A053816的公式——Select[Range[540000], Total[FromDigits /@ TakeDrop[IntegerDigits[#^2], Floor[IntegerLength[#^2]/2]]] == # &]——有问题。[code]