对“水仙花数”的研究

Study on the “Narcissistic Number”

摘要: 寻找“水仙花数”，使用实验方法，用MATLAB软件编程进行实验，得到了计算“水仙花数”的3种方法。对“水仙花数”的个数是有限的问题，用特殊值法、图解法、方程法、数学归纳法等4种方法进行了论证。对存在水仙花数的自然数位数，用MATLAB编程实验，将上界n≤60由降至n≤42。

Abstract: In order to find the “Narcissistic Number” in the natural numbers, three methods of calculating the “Narcissistic Number” are obtained by using the experimental method and MATLAB programming. The number of the “Narcissistic Number” is limited. It is demonstrated by four methods: special value method, graphic method, equation method and mathematical induction method. The digits of natural numbers contain the “Narcissistic number” that is experimented with MATLAB programming, and the upper bound is reduced from 60-digits to 42-digits.

1. 绪论

一个n位自然数，它的每个数码的n次幂之和恰好等于这个数，称这样的数为水仙花数。

当 $n=1$ 时，0至9这十个数全是水仙花数。

当 $n=2$ 时，两位自然数中没有水仙花数。事实上，假设两位自然数中有水仙花数，且设十位数为x，个位数为y，则 $x^2+y^2=10x+y$，故 $x^2-10x=-y^2+y$。由于x、y的取值范围是相同的，因此，不仿设 $y_1=x^2-10x$， $y_2=-x^2+x$。下面，用matlab编程来求解。

clear all;

clc;

clf;

x=0:1:9;

y1=x.^2-10.*x;

y2=-x.^2+x;

plot(x,y1)

hold on;

plot(x,y2)

xlabel('x');

ylabel('y');

Figure 1. The diagram of graphic method

图1. 图解法示意图

由图1可知，两曲线有两个交点 $\left(0,0\right)$， $\left(5.5,-24.75\right)$，因此，满足方程 $x^2+y^2=10x+y$ 的十位数x、个位数y不存在。故两位自然数中没有水仙花数。

当 $n=3$ 时，有四个 [1]： $153=1^3+5^3+3^3$ ； $370=3^3+7^3+0^3$ ； $371=3^3+7^3+1^3$ ； $407=4^3+0^3+7^3$。

当 $n=4$ 时，有三个： $1634=1^4+6^4+3^4+4^4$ ； $8208=8^4+2^4+0^4+8^4$ ； $9474=9^4+4^4+7^4+4^4$。

依此类推，一般地，设 $n\in N^{*}$，则存在 $M\in N$， $M=\overline{m_1{m}_2{m}_3\cdots m_n}$，使得

$m_1^n+m_2^n+m_3^n+\cdots +m_{n-1}^n+m_n^n={10}^{n-1}{m}_1+{10}^{n-2}{m}_2+{10}^{n-3}{m}_3+\cdots +10m_{n-1}+m_n$。

用matlab探索水仙花数，文献 [2] 提供了探索3位自然数中的水仙花数的matlab代码。文献 [3] 提供了一个matlab穷举法代码，自动化程度不高。但可以用于验证其他方法算出的水仙花数。这个方法，一般PC机都能验证至16位数左右。

2. 水仙花数计算方法

2.1. 计算水仙花数的方法一

对文 [2] 的方法进行改进，对输入的n位自然数分别除以10的 $n-1$ 至0次幂，余数 [4] [5] 取整得输入的自然数的各个数码，然后判断各个数码的n次幂是否等于输入的自然数，从而计算出水仙花数。因此，编制出如下matlab代码：

%求水仙花数matlab代码1.

clear all;

clc;

n=input('n=');%输入自然数的位数.

for m=10^(n-1):10^n-1%输入位数相同的自然数.

m1=rem(fix(m/10^(n-1)),10);%确定自然数的首位数码.

m2=rem(fix(m/10^(n-2)),10);%确定自然数的第二位数码，以下类推.

m3=rem(fix(m/10^(n-3)),10);

m4=rem(fix(m/10^(n-4)),10);

m5=rem(fix(m/10^(n-5)),10);

m6=rem(fix(m/10^(n-6)),10);

m7=rem(fix(m/10^(n-7)),10);

m8=rem(fix(m/10^(n-8)),10);

m9=rem(fix(m/10^(n-9)),10);

if (m1)^n+(m2)^n+(m3)^n+(m4)^n+

(m5)^n+(m6)^n+(m7)^n+(m8)^n+(m9)^n

==m%判断水仙花数

disp(m)%列出水仙花数.

end

end

这组代码，能计算出1至9位自然数中的水仙花数。如果计算机内存大，精度高 [6]，将这组代码作适当扩展，可以探索更多位数的自然数中的水仙花数。不足之处是，判断水仙花数的式子太长。

2.2. 计算水仙花数的方法二

步骤如下：第一步，输入自然数m的位数n。第二步，分离n位自然数m的各个数码。第三步，判断输入的n位自然数m是否为水仙花数。若是，则输出水仙花数。否则，返回继续搜索，直到计算出所有的水仙花数。

%求水仙花数matlab代码2.

clear all;

clc;

n=input('n=');%输入自然数m的位数n.

for m=10^(n-1):10^n-1% 输入自然数m.

x=num2str(m)-'0';%分离自然数m的各位数码.

if sum(x.^n)==m%判断一个自然数是否为水仙花数.

disp(m)%输出水仙花数.

end

end

如果计算机内存足够大，精度足够高，这个代码可以探索任意位数的自然数中的水仙花数。但是，随着位数的增加，计算用时也随着增加 [7]。

2.3. 计算水仙花数的方法三

步骤如下：第一步，输入自然数m的位数n。第二步，确定搜索自然数m的0至9这10个数码的个数的终止值(通过实验估计，取 $k=\text{fix}\left(\left(n+2\right)/3\right)$ )，以加快搜索速度。第三步，依次计算数码0至9的个数。第四步，若数码0至9的个数之和为n时，则计算各数码个数分别与其对应数码的n次幂之积的和。第五步，判断第四步所求出的和是否是n位自然数。(这样，可筛掉许多无关的数，提高了计算速度。)若是，分离其各个数码。第六步，判断这个n位自然数是否为水仙花数。若是，则输出水仙花数。否则，返回继续搜索，直到找出所有的水仙花数。

%求水仙花数matlab代码3.

clear all;

clc;

n=input('n=');%输入自然数m的位数n.

k=fix((n+2)/3);%确定搜索自然数m的0至9这10个数码的个数的终止值，加快搜索速度.

for o=0:k%搜索n位自然数m中数码0的个数.

for p=0:k%搜索n位自然数m中数码1的个数.

for q=0:k%搜索n位自然数m中数码2的个数.以下依此类推.

for r=0:k

for s=0:k

for t=0:k

for u=0:k

for v=0:k

for w=0:k

for x=0:k%搜索n位自然数m中数码9的个数.

if o+p+q+r+s+t+u+v+w+x==n%0至9这10个数码的个数的和.

m=o*0^n+p*1^n+q*2^n+r*3^n+s*4^n+t*5^n+u*6^n+v*7^n+w*8^n+x*9^n;%0至9这10个数码的n次幂的和.

if m<=10^n-1 & m>=10^(n-1)%n位自然数的范围.

y=int2str(m)-'0'; %分离自然数m的各位数码.

if sum(y.^n)==m %自然数m为水仙花数.

disp(m)%输出水仙花数.

end

end

end

end

end

end

end

end

end

end

end

end

end

如果计算机内存足够大，精度足够高，这个代码可以探索任意位数的自然数中的水仙花数，而且计算速度快。

3. 对水仙花数的个数有限问题的探究

现在有一个问题，自然数中的水仙花数的个数是有限的呢，还是无限的？就来看下面的问题：设 $n\in N^{*}$，当n为何值时，有 $n\cdot 9^n<{10}^{n-1}$ 成立 [8] [9]。

3.1. 特殊值法

事实上，对 $n\cdot 9^n<{10}^{n-1}$ 两边取自然对数，得 $\ln \left(10n\right)<n\ln \left(10/9\right)$。当 $n=60$ 时， $\ln \left(10\cdot 60\right)=6.3969$， $60\ln \left(10/9\right)=6.3216$，即有 $60\cdot 9^{60}>{10}^{60-1}$ ；当 $n=61$ 时， $\ln \left(10\cdot 61\right)=6.4135$， $61\ln \left(10/9\right)=6.4270$，即有 $61\cdot 9^{61}<{10}^{61-1}$。

3.2. 图解法

现在将离散问题变为连续函数来处理。由 $n\cdot 9^n<{10}^{n-1}$，可得 $10\cdot n<{\left(10/9\right)}^n$，将离散问题连续化，得 $10\cdot x<{\left(10/9\right)}^x$， $x\ge 1$。

解方程组 $\{\begin{array}{l}y=10x\\ y={\left(10/9\right)}^x\end{array}$ 得 $\left(60.8479,608.479\right)$。如图2所示。

Figure 2. The graph of the equations

图2. 方程组曲线图

%水仙花数图解法matlab代码.

clear all;

clf;

x=-20:0.1:70;

y1=10*x;

plot(x,y1)

hold on;

y2=exp(x*log(10/9));

plot(x,y2)

xlabel('x');

ylabel('y');

3.3. 用matlab解方程

Matlab代码为：z=fzero('10*x-exp(x*log(10/9))',[50,70]) $x=60.8479$。

因此，当自然数的位数 $n>60$ 时，就不存在水仙花数了。

3.4. 数学归纳法

将不等式 $n\cdot 9^n<{10}^{n-1}$ 等价变形为 $10\cdot n<{\left(10/9\right)}^n$。

1) 当 $n=61$ 时，有 $10\cdot 61<{\left(10/9\right)}^{61}$ ( ${\left(10/9\right)}^{61}\approx 618.3109$ )成立。2) 假设当 $n=k$ ( $k\ge 61$ )时，有 $10\cdot k<{\left(10/9\right)}^k$ 成立，那么，当 $n=k+1$ 时， ${\left(10/9\right)}^{k+1}={\left(10/9\right)}^{k}10/9>\left(10k\right)10/9=10\left(k+1\right)+\left(k-9\right)10/9>10\left(k+1\right)$ (因为 $k\ge 61$ 时，有 $\left(k-9\right)10/9>0$ )。即有 $10\left(k+1\right)<{\left(10/9\right)}^{k+1}$ 成立。由(1) (2)可知，当 $n\ge 61$，且 $n\in N^{*}$ 时，不等式 $10\cdot n<{\left(10/9\right)}^n$ 成立，亦即 $n\cdot 9^n<{10}^{n-1}$ 成立。

综上所述，当自然数的位数 $n>60$ 时，就不存在水仙花数了。因为这时，每个数码的幂指数为位数的幂之和就小于这个数了。

4. 对水仙花数的自然数位数的探究

水仙花数的个数是有限的，存在水仙花数的自然数的位数 $n\le 60$。国防科技大学刘江宁教授用计算机算出水仙花数一共有88个，其中最大的水仙花数为两个39位的自然数 [8]。

经过对系列位数较大的水仙花数各个数码的个数的分析，对存在水仙花数的自然数的位数的上界用matlab编程进行探究，可进一步降低存在水仙花数的自然数的位数的上界。

4.1. 第一次降界

对于两个39位的水仙花数，数码9都出现了7次，数码8都出现了1次，数码7都出现了4次，……表1列出了29至39位数的水仙花数的各个数码的个数。

水仙花数各个数码的个数之和等于位数n。假设n位水仙花数的数码9、8、7三个数码各有n/3个，那么，当n为何值时，有 $\frac{n}{3}\left(9^n+8^n+7^n\right)<{10}^{n-1}$，即 $9^n+8^n+7^n<\frac{3}{10n}{10}^n$ 成立？

事实上，当 $n=49$ 时， $9^{49}+8^{49}+7^{49}=5.7443e+46$， $3/10/49\cdot {10}^{49}=6.1224e+46$。即有 $9^{49}+8^{49}+7^{49}<\frac{3}{10\cdot 49}{10}^{49}$ 成立。

用matlab编程求解，z=fzero('9^x+8^x+7^x-3/10/x*10^x',[40,50]) $x=48.2515$。

因此，当自然数的位数 $n>48$ 时，就不存在水仙花数了。

Table 1. 29-39 digit narcissus number, number of digital appearance data

表1. 29-39位水仙花数各位数码出现次数数据

4.2. 第二次降界

如果假设n位水仙花数的数码9、8、7、6四个数码各有n/4个，那么，当n为何值时，有 $\frac{n}{4}\left(9^n+8^n+7^n+6^n\right)<{10}^{n-1}$，即 $9^n+8^n+7^n+6^n<\frac{2}{5n}{10}^n$ 成立？

事实上，当 $n=45$ 时， $9^{45}+8^{45}+7^{45}+6^{45}=8.7716e+42$， $2/5/45\cdot {10}^{45}=8.8889e+42$。

即有 $9^{45}+8^{45}+7^{45}+6^{45}<\frac{2}{5\cdot 45}{10}^{45}$ 成立。

用matlab编程，z=fzero('9^x+8^x+7^x+6^x-2/5/x*10^x',[40,50]) $x=44.8413$。

故当自然数的位数 $n>44$ 时，就不存在水仙花数了。

4.3. 第三次降界

如果假设n位水仙花数的数码9、8、7、6、5五个数码各有n/5个，那么，当n为何值时，有 $\frac{n}{5}\left(9^n+8^n+7^n+6^n+5^n\right)<{10}^{n-1}$，即 $9^n+8^n+7^n+6^n+5^n<\frac{1}{2n}{10}^n$ 成立？

事实上，当 $n=43$ 时， $9^{43}+8^{43}+7^{43}+6^{43}+5^{43}=1.0844e+41$， $1/2/43\cdot {10}^{43}=1.1628e+41$。

即有 $9^{43}+8^{43}+7^{43}+6^{43}+5^{43}<1/2/43\cdot {10}^{43}$ 成立。

用matlab编程，z=fzero('9^x+8^x+7^x+6^x+5^x-1/2/x*10^x',[40,50]) $x=42.1551$。

综上所述，当自然数的位数 $n>42$ 时，就不存在水仙花数了。换句话说，当自然数的位数 $n\le 42$ 时，才有可能存在水仙花数。

5. 结语

经过探讨，一是获得了用matlab计算水仙花数的3种方法；二是对存在水仙花数的自然数位数 $n\le 60$ 的问题，用4种方法进行了证明；三是将存在水仙花数的自然数的位数的上界 $n\le 60$ 经过3次降界，降至 $n\le 42$。进一步研究的问题，一是如何将存在水仙花数的自然数的位数的上界降至 $n\le 39$ ；二是水仙花数在有关学科中的应用。

致谢

感谢中国国家自然科学基金、贵州省自然科学联合基金的支持。对文中引用参考文献的作者表示衷心感谢！

基金项目

1) 国家自然科学基金项目(11661084)；2) 贵州省自然科学联合基金项目(黔科合LH字[2015]7042号)。