092. 平方数字链(Square digit chains)

将一个数的所有数字的平方相加可以得到一个新的数，不断重复这个过程，直到新的数已经出现过为止，这可以构成了一条数字链，如：
$$ \begin{align*} &44 → 32 → 13 → 10 → 1 → 1\\ &85 → 89 → 145 → 42 → 20 → 4 → 16 → 37 → 58 → 89 \end{align*} $$
可以发现，任何一个到达$1$或$89$的数字链都会陷入循环。更令人惊奇的是，从任意数开始，最终都会到达$1$或$89$。有多少个小于一千万的数最终会到达$89$？

分析：这是一道非常有趣的题目，刷过力扣的同学可能刷过202.欢乐数这道题目，题目要求我们写一个程序，判断某一个数字是否为快乐数。看一下快乐数的定义，大家就可以发现和本题的定义是一样的，只有不断计算某个数字的各位数字的平方和，如果最终到达一，那么这个数就是一个快乐数，否则就不是。在力扣的官方题解中，给出三种判断一个数是否是快乐数的方法，分别是哈希集合检测、快慢指针以及数学方法，我们重点看一下最后一种方法。

一、数学性质分析

对于平方数链，有两种情况：要不经过多次计算到达一，因为一的平方还等于一，所以链条到此时进入无限循环，我们可以认为它终止于一，也可以看作是一个不动点。另一种情况是，无论怎么计算都不会到达一，但是会进入另外一个循环：

$$ 4→16→37→58→89→145→42→20→4 $$

无论从那个数开始，只要经过若干次计算后到达这个链条中的某一个数，就会陷入到无限循环当中，题目中所说的以$89$结束的数实际上就是指计算出来的数字进入了这个循环，迟早会到达$89$。因此要判断一个平方链条是否最终会到达$89$，只需要判断计算出来的平方和是否进入这个循环就可以了。比如对于数字2，计算一次平方平方得到4，显然是循环中的数字，所以2开头的平方链条必然是以$89$结尾的。力扣的官方题解中只给出了这个观察结果，但并没有证明必然如此，会不会出现一种平方链条中的数字不断增长，即不会到达$1$也不会到达$89$的情况？下面我做一个简单的证明：

引理：对于自然数$a$，设$s(a)$表示其各位数字的平方和，则对于$a\ge100,s(a)<a$。

证：设$a$为一个$n+1$位的自然数，则据题意$n\ge2$。有：
$$ s(a)=a_n^2+a_{n-1}^2+\cdots+a_0^2 $$
因为$a_n$是一个0至9的数字，所以$a_n^2=a_na_n<10a_n$，则有：
$$ s(a)=a_n^2+a_{n-1}^2+\cdots+a_0^2<10a_n+10a_{n-1}+\cdots+10a_0 $$
又因为$a=10^na_n+10^{n-1}a_{n-1}+\cdots+a_0$，则有：
$$ a-(10a_n+10a_{n-1}+\cdots+10a_0)=(10^n-10)a_n+(10^{n-1}-10)a_{n-1}+\cdots+90a_2-9a_0 $$
因$n\ge2$，则$(10^n-10)a_n\ge90$，而$9a_0\le81$，则：
$$ (10^n-10)a_n-9a_0>0 $$
则有$a-(10a_n+10a_{n-1}+\cdots+10a_0)>0$，故：
$$ s(a)<10a_n+10a_{n-1}+\cdots+10a_0<a $$
命题得证。

因此，对于大于等于一百的自然数，计算出来的平方和必然比它自身要小，且每轮计算都会让它变得更小，因此经过多轮计算，这个数字必然跌入一百以内。下面我们只需要证明对于一百以内的自然数，$s(a)$要不最终到达$1$要不到达$89$就可以了，显然这可以通过计算机程序直接验证出来。使用图论工具，我们可以把一百以内数字的链条进行可视化，更加清晰直观一些：

在上图中每个数字都有一个箭头指向它各位数字的平方和，也就是链条中的下一个数字。我们可以非常直观的看出一百以内数字的链条分成了两种，一种是右下角的数字，最终会到达一；另一种是左上角的数字，可以看到那些以天蓝色标注的数字形成了一个环，链条中的数字一旦进入这个环就会无限循环下去。

上面我们简单证明了平方数字链条的几个数学性质，下面我们再看一下如何解决这道题。

二、快乐数的数量

使用平方数字链的性质，我们可以非常容易的验证一个数字是否是快乐数，官方题解给出的第三种算法其时间复杂度与空间复杂度分别为$O(log(n))$及$O(1)$。为了求特定范围$N$以内的快乐数或者说非快乐数的个数，我们可以一个一个数的来验证，则时间复杂度为：

$$ log(1)+log(2)+\cdots+log(N)=log(1\cdot2\cdots N)=log(N!) $$

对于规模比较大的$N$，这个时间复杂度是无法接受的。究其原因在于，在独立验证每个数是否中快乐数的过程中，涉及到非常多的重复计算，某些数字我们已经知道它是非快乐数，所以在它的链条中的数字都是非快乐数，这样可以大大缩小验证的规模。但是采用独立验证的方法，我们无法利用这个性质。

考虑到只要是三位以上的数字，它们各位数的平方和必然小于其自身，那么在经过一次平方和计算后得到数字的范围会大大缩小。对于一个$k$位数，显然是其每位数字都为9时取到最大，因此其平方和最大为$9^2k=81k$。比如当$k=7$时，最大的平方和为$567$。显然，对于一个$k$位数，经过一轮计算，如果计算出来的数字是一个快乐数，那个这个数字就一定是快乐数，反之就是一个非快乐数。

现在我们可以反过来思考，如果一个数字是快乐数，那么有那些数字在经过一轮计算后会得到这个数。或者更一般地，我们设$f(n,k)$表示$k$位以内的数字中其平方和等于$n$的数字的个数，我们应该如何求解$f(n,k)$。设$k$位数的每位数字为$d_i(1\le i\le k)$，则有：

$$ d_1^2+d_2^2+\cdots+d_k^2=n\Rightarrow d_1^2+d_2^2+\cdots+d_{k-1}^2=n-d_k^2 $$

为了求上式右侧满足要求的个数，我们可以把$d_k=0,1,2\cdots,9$依次代入，然后看$d_1$至$d_{k-1}$的平方和等于$n,n-1^2,n-2^2,\cdots,n-9^2$的个数，而上式右侧实际上可以看作是把$n-d_k^2$表示成为$k-1$位数的平方和的方法数，即为$f(n-d_k^2,k-1)，d_k=0,1,2\cdots,9$，则有：

$$ f(n,k)=\sum_{i=0}^9f(n-i^2,k-1) $$

根据以上递推式，我们可以用动态规划的方法来求解$f(n,k)$。边界条件有三个：

如果$n<0$，一个负数显示无法表示为平方和，则此时$f(n,k)=0$；
如果$n>0$且$k=0$，一个正数无法表示成零个数的平方和，则此时$f(n,k)=0$；
如果$n=k=0$，即用零个数来表示零，此时$f(n,k)=1$。

因此完整的状态转移方程为：

$$ f(n,k)=\begin{cases} 0 & \text{ if } n<0 \\ 0 & \text{ if } n>0\ \&\ k=0 \\ 1 & \text{ if } n=k=0 \\ \sum_{i=0}^9f(n-i^2,k-1) & \text{ else } \end{cases} $$

三、代码实现

综合前面两部分的分析，为了计算$10^7$内的非快乐数的数量，我们只需要计算出该范围内的快乐数数量，再从总数里减去就可以了。对于$10^7$内的数字，其各位数字平方和最大为$567$，因此我们只需要找到这个范围内的快乐数即可，因为数据规模较小，我们只需要逐一判断即可。在找到这个范围内的快乐数后，对于每一个快乐数$n$，我们计算$f(n,7)$，再将结果加总起来即为$10^7$内的快乐数的数量。代码如下：

# time cost = 1.95 ms ± 12.1 µs

from functools import lru_cache
from math import log10

def is_happy(n):

    cycle_members = {4, 16, 37, 58, 89, 145, 42, 20}

    def get_next(number):
        total_sum = 0
        while number > 0:
            number, digit = divmod(number, 10)
            total_sum += digit ** 2
        return total_sum

    while n != 1 and n not in cycle_members:
        n = get_next(n)

    return n == 1

@lru_cache(maxsize=None)
def f(n,k):
    if n == k == 0:
        return 1
    elif n < 0:
        return 0
    elif n > 0 and k ==0:
        return 0
    else:
        return sum([f(n-i**2,k-1) for i in range(10)])

def main(N=10**7):
    power = int(log10(N))
    u = 9**2*power
    happy_numbers = [i for i in range(1,u+1) if is_happy(i)]
    ways = sum([f(x,power) for x in happy_numbers])
    return (N - ways - 1)

使用以上算法，我们甚至可以在一秒时间内计算出$10^{100}$内的快乐数的数量。

四、相关延伸

快乐数是通过计算十进制下的各位数字的平方和来判断的，这一概念可以进行类比推广。一般地，对于$b$进制下数字，计算其每位数的$p$次方之和，如果这个和与这个数字相等，我们称这种数为超完全数字不变数，有时候也叫水仙花数或者自恋数。可以证明，对于数$n$在$b$进制下各位数字的$p$次方之和$S_{p,b}(n)$，当${\displaystyle n\leq (p-2)^{p}+p(b-1)^{p}}$时，$S_{p,b}(n)$要不会变成一个不动点，要不会进入一个循环，具体的证明大家可以参见维基百科。对于$p=2,b=10$的特殊情况，也就是我们在前文中得到的结果。对于十进制下的其它次方，一个简要的结果如下：

微信截图_20210302151548

五、参考

本站所有文章如非特别声明皆为原创，转载请与作者本人联系并注明出处(原文链接)。如对文章有任何疑问或需指出错误，可以此页面提出问题，我会尽量及时回复。

←返回主页