关于幂等差型水仙花数的研究

韦佳玉，杨健，蒋剑军

（铜陵学院数学与计算机学院，铜陵 244061）

0 引言

水仙花数是一类广为人知的正整数。最初的水仙花数是指一个3位数，它各位数字的立方和等于该数自身［1］。例如，153就是一个水仙花数，因为153=13+53+33。众所周知，3位水仙花数有4个：153、370、371和407。

水仙花数是如此的有趣，以致在众多的科普学者中口口相传，也成为编程爱好者推崇的经典案例。比如，谈祥柏在文献［2］中对水仙花数、黑洞与振荡方面的知识做了科普；马丽燕等在文献［3］中以水仙花数为例介绍了Python的程序设计方法；史婉玉等则在文献［4］中介绍了水仙花数在Java程序设计中的应用。有些学者则从水仙花数挖掘编程的内涵［5］或编程课程对学习习惯的培养［6］。水仙花数甚至渗透到了高中信息技术的学习中［7］。

3位水仙花数是程序设计中的优秀案例。不仅如此，它也是理论探索的一个源泉。林宣治［8］将水仙花数推广到高位的情况。所谓高位水仙花数，是指一个n位正整数，它各位数字的n次方和等于该数自身。卫洪春［9］则设计了基于C语言动态数组的快速求解算法计算出了所有的高位水仙花数。

近期，杨健在文献［10］中将水仙花数从高位情形推广到指数位移情形。所谓位移水仙花数，是指一个n位正整数，它能表示为它各位数字的n+d次幂之和，即：

其中d(≥-1)为整数，称为位移。

本文将经典水仙花数［1］、高位水仙花数［8］及位移水仙花数［10］统称为齐次水仙花数。

继3位水仙花数之后，寻找齐次水仙花数成为编程进阶的优秀案例，有着独特的意义。一方面是计算机程序设计教学中典型的程序设计范例，因为寻找水仙花数包含着基础编程所必需掌握的选择结构与循环结构；另一方面又是整数理论中一类具有独特魅力的整数，引出了齐次水仙花数型的不定方程，即形如

的不定方程，其中xi(i=1,…,n)都是介于0到9的整数，且x1≠0，k为待定次数。文献［10］提出了关于齐次水仙花数和不定方程（1）的两个公开问题。

本文将齐次水仙花数的概念做进一步的推广，提出幂等差型水仙花数。首先给出幂等差型水仙花数的定义，然后讨论了幂等差型水仙花数的性质，最后设计了基于Matlab的快速求解算法。

1 幂等差型水仙花数

1.1 定义

设b和d为两个整数，且b≥0。若一个n位正整数，能表示为它各位数字的如下形式的幂之和：

则称该n位正整数为基数为b公差为d的幂等差型水仙花数。b称为指数序列的基数；d称为指数序列的公差。

式（2）中 当d≤0时，令c=-d≥0,a=b+(n-1)d≥0，则（2）式变为：

所以，式（2）中当d≥0时称为幂增型水仙花数，当d≤0式称为幂减型水仙花数。式（3）正是幂减型水仙花数的直观形式。

当d=0时的幂等差水仙花数正是齐次水仙花数，所以文献［10］研究的位移水仙花数是本文的一个特例。

通过计算知，598=51+92+83是基为b=1公差为d=1的幂增数；332=35+34+23是基为b=3公差为d=-1的幂减数；153=13+53+33是b=3的齐次数。

因为位上数字可能有1、0或重复等情形，整数的幂等差表示可能不唯一，比如

所以1676既是幂增数又是幂减数，原因在于位上数字出现了1和二重数字6。

1.2 性质

因为幂等差水仙花数有两个参数：数b和d，所以具有比齐次水仙花数更丰富和更复杂从而更难以揭秘的内蕴性质。

性质1有限性条件。设b和d为两个给定整数且b≥0。则当|d|≤1时，基为b公差为d的幂等差型水仙花数是有限的。

证明设是基数为b公差为d的幂等差数，则式（2）成立。分三种情形证明在题设条件下幂等差数的有限性。

情形1：d=0。已由文献［10］证明。

情形2：d>0。则式（2）右端可放大为：

由式（4）进而得：

式（5）首尾两端取以10为底的对数，得：

式（6）表明，当d≤1，即d=1时，有：

情形3：d<0。则令c=-d≥0,a=b-(n+1)c≥0，由式（2）变为式（3）；对式（3）右端进行如式（4）一样的放大，得到：

进而得到c=1及与式（7）一致的表达式

注：性质1表明，对于给定的b(≥1)，||d≤1是“基为b公差为d的幂等差水仙花数是有限的”的充分条件。由于式（4）右端放得过大，|d|≤1不可能成为有限性的必要条件。因为公差|d|>1的幂等差数是存在的，比如，463=41+63+35是公差为d=2的幂增数，所以下面两个问题是自然而然的：

（1）幂等差数的公差d是否随着位数n的变化而有界？

（2）对于任意给定的两个整数b(≥1)和d，基为b公差为d的幂等差数是否有限？

性质2n与b和d的关系性质。当||d≤1时，幂等差型水仙花数的位数n与基b和公差d之间满足如下关系式：

d值0±1关系式n-lg n-1 lg 9

证明（10）式已由文献［10］给出，（11）式由式（7）给出。

应用性质2的关系式，由Matlab编程计算给定d,b对应的n值的范围如表1和表2所示：

表1 d=0时n随b变化的下界和上界

表2 d=±1时n随b变化的上界

注：应用式（10）和式（11）得到的n的上界和下界，并非n的上确界和下确界。

性质3位异性。假设是 基 数为b公差为d的幂等差型水仙花数，则

（1）位上的数字必有不同于0或1的数字；

证明：反证法先证明（1）。

设n位由0和1构成的正整数是基数为b公差为d的幂等差数，则

又易知，当且仅当n=1,时，

再反证法证明（2）。

设n位正整数是基数为b公差为d的幂等差数。则m是1到9的自然数,由上述1的证明可知，m一定不为1。

于是，有：

又

结合上式，得到：

故

可知当且仅当b=d=1,m=10时，等式成立。与m是2到9的自然数矛盾。故至少有两个位上的数字不同。

2 幂等差型水仙花数的快速求解算法

计算幂等差型水仙花数，计算量超大是其典型的特点。文献［9］设计了高位水仙花数的C语言快速求解算法。本文设计基于Matlab的快速计算泛水仙花数的算法。在设计基于Matlab实现的算法时，须充分考虑Matlab的计算特点。众所周知，Matlab计算特点是数据以矩阵为单元，循环效率相对较低。所以在设计算法时必须将输入、输出和存储数据都矩阵化，降低对循环的依赖，以提高效率。

2.1 算法描述

算法的核心思想就是基于Matlab以矩阵为计算单元，矩阵化中间计算数据。

（1）将位数为n的所有整数矩阵化，放在变量allnumbers中，即:

在Matlab中这是一个行向量。

所有n位数有k=9*10n-1个。

（2）用Matlab命令套装num2str、str2num及reshape提取allnumbers中所有n位整数各个位上的数字，构成一个n×k矩阵，放在变量bitsMa⁃trix中，即:

上述三个语句，第一个语句得到的a是一个字符串行向量，长度为k(n+2)-2；第二个语句是将第一个语句的结果转置后转化为整数，所得结果b是一个长度为nk的列向量；第三个语句是将b转换为n×k矩阵，该矩阵的第j列就是提取出来的第j个n位整数各位上的数字。

（3）将指数序列base+d,…,base+nd转化为n×k矩阵，放在变量e中，实现语句如下

通过上述语句，即得n×k指数矩阵e。

（4）计算矩阵bitsMatrix的e次幂，并按列求和，结果放在变量esum中，实现语句为：

（5）判断向量allnumbers与esum中对应位置的数是否相等：若相等，是幂等差数；若不等，则不是幂等差水仙花数。

2.2 计算结果

由于9位数以上的数据利用Matlab运算时间过长，我们在这里只计算了9位数以内的幂增型泛水仙花数。

9位数以内的基数为b的幂增型泛水仙花数,且显然位数为1时只有1是泛水仙花数，故不再计算位数为1的情况。

给定位数n，计算8位数以内的所有幂增型泛水仙花数需要7141.182828秒即约为119分钟的时间。表3列出了108以内的幂等差水仙花数。

表3 108以内的幂等差数

3 结语

本文将林宣治、杨健等研究的齐次水仙花数推广为更一般的幂等差型水仙花数，使得齐次水仙花数成为我们研究范畴的一个特例，并对幂等差型水仙花数的性质进行了研究。

在对幂等差型水仙花数研究的过程中，不禁产生了一些疑惑，其中最大的疑惑就是：

问题1是否存在任意公差的幂等差型水仙花数？

显然，无论计算机技术发展到什么程度都无法从计算的角度解决上述问题。于是，我们把上述问题等价地转化为如下数学问题：

问题2下述关于(n;b,d,x1,x2,…,xn)的不定方程

的解是否是有限的？其中0≤xk≤9,k=1,2,…,n，且xn≠0。

期待数学工作者给出问题2的解答。