计算机中的数学【水仙花数】求解自然数中所有的水仙花数-阿里云开发者社区

计算机中的数学【水仙花数】求解自然数中所有的水仙花数

2018-05-06 2014

版权

本文内容由阿里云实名注册用户自发贡献，版权归原作者所有，阿里云开发者社区不拥有其著作权，亦不承担相应法律责任。具体规则请查看《阿里云开发者社区用户服务协议》和《阿里云开发者社区知识产权保护指引》。如果您发现本社区中有涉嫌抄袭的内容，填写侵权投诉表单进行举报，一经查实，本社区将立刻删除涉嫌侵权内容。

简介： 在数论中，水仙花数（Narcissistic number），也被称为超完全数字不变数（pluperfect digital invariant, PPDI）、自恋数、自幂数、阿姆斯壮数或阿姆斯特朗数（Armstrong number），用来描述一个N位非负整数，其各位数字的N次方和等于该数本身。

在数论中，水仙花数（Narcissistic number），也被称为超完全数字不变数（pluperfect digital invariant, PPDI）、自恋数、自幂数、阿姆斯壮数或阿姆斯特朗数（Armstrong number），用来描述一个N位非负整数，其各位数字的N次方和等于该数本身。

水仙花数只是自幂数的一种，严格来说3位数的3次幂数才称为水仙花数。
附：其他位数的自幂数名字
一位自幂数：独身数
两位自幂数：没有
三位自幂数：水仙花数
四位自幂数：四叶玫瑰数
五位自幂数：五角星数
六位自幂数：六合数
七位自幂数：北斗七星数
八位自幂数：八仙数
九位自幂数：九九重阳数
十位自幂数：十全十美数

可以证明：

当 n > 60 的时候，有

也就是说，n 位水仙花数最小的数字是10^(n-1) , 例如，3位水仙花数最小是10^2 = 100, 这个是个，n 位最小的数都大于各个位上的数字的 n 幂次和最大值： n * 9^n 。

最大的水仙花数有39位。十进制自然数中的所有水仙花数共有88个。

image.png

使用 Kotlin 编程来计算自然数中所有的水仙花数。

完整代码如下：

package com.easykotlin.lectures.kfp.math

import java.math.BigInteger

class NarcissisticNumbers {
    val TEN = BigInteger.TEN

    /**
     *
     * 获取一个 n 位数 N 的各个位上的数字，放到一个 List 中：
     * N=153, digits=[1,5,3]
     */
    fun digits(N: BigInteger, n: Int): ArrayList<Int> {
        var M = N
        var digits = arrayListOf<Int>()

        (1..n).forEach {
            val remainder = M.mod(TEN)
            digits.add(remainder.intValueExact())
            M = M.divide(TEN)
        }
        digits.reverse()
        return digits
    }

    /**
     *  数字 0-9 的 n次幂
     *  0, 1, 2^n, 3^n,...,9^n
     */
    fun zero2NinePower(n: Int): List<BigInteger> {
        var result = arrayListOf<BigInteger>()
        result.add(BigInteger.ZERO)
        result.add(BigInteger.ONE)
        (2..9).forEach {
            result.add(BigInteger.valueOf(it.toLong()).pow(n))
        }
        return result
    }

    fun checkOut(n: Int) {
        var N_ = BigInteger.ZERO
        val zero2NinePowers = zero2NinePower(n)
        // d0 表示0出现次数，d1 表示1出现次数。例如： 数字153，有 d0=1,d1=1,d2=0,d3=1,d4=0,d5=1,...
        for (d0 in 0L..(n - 1)) {
            for (d1 in 0L..(n - d0)) {
                for (d2 in 0L..(n - d0 - d1)) {
                    for (d3 in 0L..(n - d0 - d1 - d2)) {
                        for (d4 in 0L..(n - d0 - d1 - d2 - d3)) {
                            for (d5 in 0L..(n - d0 - d1 - d2 - d3 - d4)) {
                                for (d6 in 0L..(n - d0 - d1 - d2 - d3 - d4 - d5)) {
                                    for (d7 in 0L..(n - d0 - d1 - d2 - d3 - d4 - d5 - d6)) {
                                        for (d8 in 0L..(n - d0 - d1 - d2 - d3 - d4 - d5 - d6 - d7)) {
                                            for (d9 in 0L..(n - d0 - d1 - d2 - d3 - d4 - d5 - d6 - d7 - d8)) {
                                                N_ = zero2NinePowers[0].multiply(BigInteger.valueOf(d0)) +
                                                        zero2NinePowers[1].multiply(BigInteger.valueOf(d1)) +
                                                        zero2NinePowers[2].multiply(BigInteger.valueOf(d2)) +
                                                        zero2NinePowers[3].multiply(BigInteger.valueOf(d3)) +
                                                        zero2NinePowers[4].multiply(BigInteger.valueOf(d4)) +
                                                        zero2NinePowers[5].multiply(BigInteger.valueOf(d5)) +
                                                        zero2NinePowers[6].multiply(BigInteger.valueOf(d6)) +
                                                        zero2NinePowers[7].multiply(BigInteger.valueOf(d7)) +
                                                        zero2NinePowers[8].multiply(BigInteger.valueOf(d8)) +
                                                        zero2NinePowers[9].multiply(BigInteger.valueOf(d9))
                                                doCheck(N_, d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, n)
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }

    private fun doCheck(N_: BigInteger, d0: Long, d1: Long, d2: Long, d3: Long, d4: Long, d5: Long, d6: Long, d7: Long, d8: Long, d9: Long, n: Int) {
        // 数字 N_ 每位上的数字，例如： N_ = 153，  digitCounts = [1,5,3]
        val digitCounts = digits(N_, n)
        // 统计 digitCounts 中 0-9 分别出现的次数
        val d_ = arrayListOf(0L, 0, 0, 0, 0, 0, 0, 0, 0, 0)
        digitCounts.forEach {
            when (it) {
                0 -> d_[0]++
                1 -> d_[1]++
                2 -> d_[2]++
                3 -> d_[3]++
                4 -> d_[4]++
                5 -> d_[5]++
                6 -> d_[6]++
                7 -> d_[7]++
                8 -> d_[8]++
                9 -> d_[9]++
            }
        }
        var sum = 0
        d_.forEach { sum += it.toInt() }
        if (d0 == d_[0] &&
                d1 == d_[1] &&
                d2 == d_[2] &&
                d3 == d_[3] &&
                d4 == d_[4] &&
                d5 == d_[5] &&
                d6 == d_[6] &&
                d7 == d_[7] &&
                d8 == d_[8] &&
                d9 == d_[9] &&
                sum == n && validateBitWidth(N_, n)) { // 完全数字不变数
            println("${N_}")
//            println("d_ = ${d_}")
//            println("digitCounts = ${digitCounts}")
        }
    }

    /** 过滤掉，N 首位是 0 的情况：
     * 类似 : 5 位数：
    这样的数称为完全数字不变数（perfect digital invariant）
    N_ = 4151
    d_ = [1, 2, 0, 0, 1, 1, 0, 0, 0, 0]
    digitCounts = [0, 4, 1, 5, 1]
    5 位数：
    N_ = 4150
    d_ = [2, 1, 0, 0, 1, 1, 0, 0, 0, 0]
    digitCounts = [0, 4, 1, 5, 0]
    5 位数：
    N_ = 1
    d_ = [4, 1, 0, 0, 0, 0, 0, 0, 0, 0]
    digitCounts = [0, 0, 0, 0, 1]
     */
    fun validateBitWidth(N_: BigInteger, n: Int): Boolean {
        var count = 0
        var N = N_
        while (N > BigInteger.ZERO) {
            N = N.divide(BigInteger.TEN)
            count++
        }
//        println("count=$count")
        return n == count
    }

}

执行入口函数

    @Test
    fun test_checkOut() {
        (1..39).forEach {
            println("-----------------------------")
            println("$it 位水仙花数：")
            val s = System.currentTimeMillis()
            NarcissisticNumbers.checkOut(it)
            val t = System.currentTimeMillis()
            println("时间：${t - s} ms")
        }
    }

运行结果：

-----------------------------
1 位水仙花数：
9
8
7
6
5
4
3
2
1
时间：42 ms
-----------------------------
2 位水仙花数：
时间：8 ms
-----------------------------
3 位水仙花数：
371
153
407
370
时间：21 ms
-----------------------------
4 位水仙花数：
9474
1634
8208
时间：45 ms
-----------------------------
5 位水仙花数：
54748
92727
93084
时间：35 ms
-----------------------------
6 位水仙花数：
548834
时间：190 ms
-----------------------------
7 位水仙花数：
9926315
1741725
4210818
9800817
时间：142 ms
-----------------------------
8 位水仙花数：
88593477
24678051
24678050
时间：415 ms
-----------------------------
9 位水仙花数：
534494836
472335975
912985153
146511208
时间：468 ms
-----------------------------
10 位水仙花数：
4679307774
时间：1428 ms
-----------------------------
11 位水仙花数：
82693916578
44708635679
94204591914
32164049651
49388550606
42678290603
40028394225
32164049650
时间：1374 ms
-----------------------------
12 位水仙花数：
时间：2187 ms
-----------------------------
13 位水仙花数：
时间：4884 ms
-----------------------------
14 位水仙花数：
28116440335967
时间：8053 ms
-----------------------------
15 位水仙花数：
时间：14886 ms
-----------------------------
16 位水仙花数：
4338281769391371
4338281769391370
时间：24496 ms
-----------------------------
17 位水仙花数：
35641594208964132
21897142587612075
35875699062250035
时间：40413 ms
-----------------------------
18 位水仙花数：
时间：70124 ms
-----------------------------
19 位水仙花数：
4498128791164624869
4929273885928088826
3289582984443187032
1517841543307505039
时间：93881 ms
-----------------------------
20 位水仙花数：
63105425988599693916
时间：159871 ms
-----------------------------
21 位水仙花数：
128468643043731391252
449177399146038697307
时间：244321 ms
-----------------------------
22 位水仙花数：
时间：395756 ms
-----------------------------
23 位水仙花数：
21887696841122916288858
28361281321319229463398
27879694893054074471405
35452590104031691935943
27907865009977052567814
时间：612641 ms
-----------------------------
24 位水仙花数：
188451485447897896036875
239313664430041569350093
174088005938065293023722
时间：888945 ms
-----------------------------
25 位水仙花数：
...

使用一台普通的 PC 机器（单机、单线程）：

可以看出——
前15位水仙花数，在 10 s 时间量级；
21位水仙花数，时间 4 min 。
22位数字中没有水仙花数。花费 5min。
23位水仙花数，时间 10 min 。
24位水仙花数，时间 15 min 。
......
后面的位数越大，时间将会翻倍。不过，终归会在有限的天数内完成计算。

当然，现代超大规模、并行计算机算起来会快很多。
上面的算法也有进一步优化的空间。

算法代码中的函数说明如下：

zero2NinePower() 函数：

    @Test
    fun test_zero2NinePower() {
        println(NarcissisticNumbers.zero2NinePower(3))
        println(NarcissisticNumbers.zero2NinePower(4))
        println(NarcissisticNumbers.zero2NinePower(21))
    }

输出：

[0, 1, 8, 27, 64, 125, 216, 343, 512, 729]
[0, 1, 16, 81, 256, 625, 1296, 2401, 4096, 6561]
[0, 1, 2097152, 10460353203, 4398046511104, 476837158203125, 21936950640377856, 558545864083284007, 9223372036854775808, 109418989131512359209]

digits() 函数

    @Test
    fun test_digits() {
        println(NarcissisticNumbers.digits(BigInteger.valueOf(153), 3))
        println(NarcissisticNumbers.digits(BigInteger.valueOf(1634), 4))
        println(NarcissisticNumbers.digits(BigInteger.valueOf(4150), 4))
    }

输出：