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\)就可以了,显然这可以通过计算机程序直接验证出来。使用图论工具,我们可以把一百以内数字的链条进行可视化,更加清晰直观一些:

2

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

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

二、快乐数的数量

使用平方数字链的性质,我们可以非常容易的验证一个数字是否是快乐数,官方题解给出的第三种算法其时间复杂度与空间复杂度分别为\(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)\)。边界条件有三个:

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

$$ 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

五、参考

  1. Todd, Elizabeth, "Happy numbers" (2018). Creative Components . 109.
  2. Gilmer, J. (2011). On the density of happy numbers. arXiv preprint arXiv:1110.3836.
  3. Wikipedia contributors. (2021, February 1). Happy number. In Wikipedia
  4. Wikipedia contributors. (2020, November 2). Perfect digital invariant. In Wikipedia