韋佳玉，楊健，蔣劍軍

（銅陵學(xué)院數(shù)學(xué)與計(jì)算機(jī)學(xué)院，銅陵 244061）

0 引言

水仙花數(shù)是一類廣為人知的正整數(shù)。最初的水仙花數(shù)是指一個(gè)3位數(shù)，它各位數(shù)字的立方和等于該數(shù)自身［1］。例如，153就是一個(gè)水仙花數(shù)，因?yàn)?53=13+53+33。眾所周知，3位水仙花數(shù)有4個(gè)：153、370、371和407。

水仙花數(shù)是如此的有趣，以致在眾多的科普學(xué)者中口口相傳，也成為編程愛好者推崇的經(jīng)典案例。比如，談祥柏在文獻(xiàn)［2］中對水仙花數(shù)、黑洞與振蕩方面的知識做了科普；馬麗燕等在文獻(xiàn)［3］中以水仙花數(shù)為例介紹了Python的程序設(shè)計(jì)方法；史婉玉等則在文獻(xiàn)［4］中介紹了水仙花數(shù)在Java程序設(shè)計(jì)中的應(yīng)用。有些學(xué)者則從水仙花數(shù)挖掘編程的內(nèi)涵［5］或編程課程對學(xué)習(xí)習(xí)慣的培養(yǎng)［6］。水仙花數(shù)甚至滲透到了高中信息技術(shù)的學(xué)習(xí)中［7］。

3位水仙花數(shù)是程序設(shè)計(jì)中的優(yōu)秀案例。不僅如此，它也是理論探索的一個(gè)源泉。林宣治［8］將水仙花數(shù)推廣到高位的情況。所謂高位水仙花數(shù)，是指一個(gè)n位正整數(shù)，它各位數(shù)字的n次方和等于該數(shù)自身。衛(wèi)洪春［9］則設(shè)計(jì)了基于C語言動態(tài)數(shù)組的快速求解算法計(jì)算出了所有的高位水仙花數(shù)。

近期，楊健在文獻(xiàn)［10］中將水仙花數(shù)從高位情形推廣到指數(shù)位移情形。所謂位移水仙花數(shù)，是指一個(gè)n位正整數(shù)，它能表示為它各位數(shù)字的n+d次冪之和，即：

其中d(≥-1)為整數(shù)，稱為位移。

本文將經(jīng)典水仙花數(shù)［1］、高位水仙花數(shù)［8］及位移水仙花數(shù)［10］統(tǒng)稱為齊次水仙花數(shù)。

繼3位水仙花數(shù)之后，尋找齊次水仙花數(shù)成為編程進(jìn)階的優(yōu)秀案例，有著獨(dú)特的意義。一方面是計(jì)算機(jī)程序設(shè)計(jì)教學(xué)中典型的程序設(shè)計(jì)范例，因?yàn)閷ふ宜苫〝?shù)包含著基礎(chǔ)編程所必需掌握的選擇結(jié)構(gòu)與循環(huán)結(jié)構(gòu)；另一方面又是整數(shù)理論中一類具有獨(dú)特魅力的整數(shù)，引出了齊次水仙花數(shù)型的不定方程，即形如

的不定方程，其中xi(i=1,…,n)都是介于0到9的整數(shù)，且x1≠0，k為待定次數(shù)。文獻(xiàn)［10］提出了關(guān)于齊次水仙花數(shù)和不定方程（1）的兩個(gè)公開問題。

本文將齊次水仙花數(shù)的概念做進(jìn)一步的推廣，提出冪等差型水仙花數(shù)。首先給出冪等差型水仙花數(shù)的定義，然后討論了冪等差型水仙花數(shù)的性質(zhì)，最后設(shè)計(jì)了基于Matlab的快速求解算法。

1 冪等差型水仙花數(shù)

1.1 定義

設(shè)b和d為兩個(gè)整數(shù)，且b≥0。若一個(gè)n位正整數(shù)，能表示為它各位數(shù)字的如下形式的冪之和：

則稱該n位正整數(shù)為基數(shù)為b公差為d的冪等差型水仙花數(shù)。b稱為指數(shù)序列的基數(shù)；d稱為指數(shù)序列的公差。

式（2）中 當(dāng)d≤0時(shí)，令c=-d≥0,a=b+(n-1)d≥0，則（2）式變?yōu)椋?/p>

所以，式（2）中當(dāng)d≥0時(shí)稱為冪增型水仙花數(shù)，當(dāng)d≤0式稱為冪減型水仙花數(shù)。式（3）正是冪減型水仙花數(shù)的直觀形式。

當(dāng)d=0時(shí)的冪等差水仙花數(shù)正是齊次水仙花數(shù)，所以文獻(xiàn)［10］研究的位移水仙花數(shù)是本文的一個(gè)特例。

通過計(jì)算知，598=51+92+83是基為b=1公差為d=1的冪增數(shù)；332=35+34+23是基為b=3公差為d=-1的冪減數(shù)；153=13+53+33是b=3的齊次數(shù)。

因?yàn)槲簧蠑?shù)字可能有1、0或重復(fù)等情形，整數(shù)的冪等差表示可能不唯一，比如

所以1676既是冪增數(shù)又是冪減數(shù)，原因在于位上數(shù)字出現(xiàn)了1和二重?cái)?shù)字6。

1.2 性質(zhì)

因?yàn)閮绲炔钏苫〝?shù)有兩個(gè)參數(shù)：數(shù)b和d，所以具有比齊次水仙花數(shù)更豐富和更復(fù)雜從而更難以揭秘的內(nèi)蘊(yùn)性質(zhì)。

性質(zhì)1有限性條件。設(shè)b和d為兩個(gè)給定整數(shù)且b≥0。則當(dāng)|d|≤1時(shí)，基為b公差為d的冪等差型水仙花數(shù)是有限的。

證明設(shè)是基數(shù)為b公差為d的冪等差數(shù)，則式（2）成立。分三種情形證明在題設(shè)條件下冪等差數(shù)的有限性。

情形1：d=0。已由文獻(xiàn)［10］證明。

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

由式（4）進(jìn)而得：

式（5）首尾兩端取以10為底的對數(shù)，得：

式（6）表明，當(dāng)d≤1，即d=1時(shí)，有：

情形3：d<0。則令c=-d≥0,a=b-(n+1)c≥0，由式（2）變?yōu)槭剑?）；對式（3）右端進(jìn)行如式（4）一樣的放大，得到：

進(jìn)而得到c=1及與式（7）一致的表達(dá)式

注：性質(zhì)1表明，對于給定的b(≥1)，||d≤1是“基為b公差為d的冪等差水仙花數(shù)是有限的”的充分條件。由于式（4）右端放得過大，|d|≤1不可能成為有限性的必要條件。因?yàn)楣顋d|>1的冪等差數(shù)是存在的，比如，463=41+63+35是公差為d=2的冪增數(shù)，所以下面兩個(gè)問題是自然而然的：

（1）冪等差數(shù)的公差d是否隨著位數(shù)n的變化而有界？

（2）對于任意給定的兩個(gè)整數(shù)b(≥1)和d，基為b公差為d的冪等差數(shù)是否有限？

性質(zhì)2n與b和d的關(guān)系性質(zhì)。當(dāng)||d≤1時(shí)，冪等差型水仙花數(shù)的位數(shù)n與基b和公差d之間滿足如下關(guān)系式：

d值0±1關(guān)系式n-lg n-1 lg 9

證明（10）式已由文獻(xiàn)［10］給出，（11）式由式（7）給出。

應(yīng)用性質(zhì)2的關(guān)系式，由Matlab編程計(jì)算給定d,b對應(yīng)的n值的范圍如表1和表2所示：

表1 d=0時(shí)n隨b變化的下界和上界

表2 d=±1時(shí)n隨b變化的上界

注：應(yīng)用式（10）和式（11）得到的n的上界和下界，并非n的上確界和下確界。

性質(zhì)3位異性。假設(shè)是 基 數(shù)為b公差為d的冪等差型水仙花數(shù)，則

（1）位上的數(shù)字必有不同于0或1的數(shù)字；

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

設(shè)n位由0和1構(gòu)成的正整數(shù)是基數(shù)為b公差為d的冪等差數(shù)，則

又易知，當(dāng)且僅當(dāng)n=1,時(shí)，

再反證法證明（2）。

設(shè)n位正整數(shù)是基數(shù)為b公差為d的冪等差數(shù)。則m是1到9的自然數(shù),由上述1的證明可知，m一定不為1。

于是，有：

又

結(jié)合上式，得到：

故

可知當(dāng)且僅當(dāng)b=d=1,m=10時(shí)，等式成立。與m是2到9的自然數(shù)矛盾。故至少有兩個(gè)位上的數(shù)字不同。

2 冪等差型水仙花數(shù)的快速求解算法

計(jì)算冪等差型水仙花數(shù)，計(jì)算量超大是其典型的特點(diǎn)。文獻(xiàn)［9］設(shè)計(jì)了高位水仙花數(shù)的C語言快速求解算法。本文設(shè)計(jì)基于Matlab的快速計(jì)算泛水仙花數(shù)的算法。在設(shè)計(jì)基于Matlab實(shí)現(xiàn)的算法時(shí)，須充分考慮Matlab的計(jì)算特點(diǎn)。眾所周知，Matlab計(jì)算特點(diǎn)是數(shù)據(jù)以矩陣為單元，循環(huán)效率相對較低。所以在設(shè)計(jì)算法時(shí)必須將輸入、輸出和存儲數(shù)據(jù)都矩陣化，降低對循環(huán)的依賴，以提高效率。

2.1 算法描述

算法的核心思想就是基于Matlab以矩陣為計(jì)算單元，矩陣化中間計(jì)算數(shù)據(jù)。

（1）將位數(shù)為n的所有整數(shù)矩陣化，放在變量allnumbers中，即:

在Matlab中這是一個(gè)行向量。

所有n位數(shù)有k=9*10n-1個(gè)。

（2）用Matlab命令套裝num2str、str2num及reshape提取allnumbers中所有n位整數(shù)各個(gè)位上的數(shù)字，構(gòu)成一個(gè)n×k矩陣，放在變量bitsMa?trix中，即:

上述三個(gè)語句，第一個(gè)語句得到的a是一個(gè)字符串行向量，長度為k(n+2)-2；第二個(gè)語句是將第一個(gè)語句的結(jié)果轉(zhuǎn)置后轉(zhuǎn)化為整數(shù)，所得結(jié)果b是一個(gè)長度為nk的列向量；第三個(gè)語句是將b轉(zhuǎn)換為n×k矩陣，該矩陣的第j列就是提取出來的第j個(gè)n位整數(shù)各位上的數(shù)字。

（3）將指數(shù)序列base+d,…,base+nd轉(zhuǎn)化為n×k矩陣，放在變量e中，實(shí)現(xiàn)語句如下

通過上述語句，即得n×k指數(shù)矩陣e。

（4）計(jì)算矩陣bitsMatrix的e次冪，并按列求和，結(jié)果放在變量esum中，實(shí)現(xiàn)語句為：

（5）判斷向量allnumbers與esum中對應(yīng)位置的數(shù)是否相等：若相等，是冪等差數(shù)；若不等，則不是冪等差水仙花數(shù)。

2.2 計(jì)算結(jié)果

由于9位數(shù)以上的數(shù)據(jù)利用Matlab運(yùn)算時(shí)間過長，我們在這里只計(jì)算了9位數(shù)以內(nèi)的冪增型泛水仙花數(shù)。

9位數(shù)以內(nèi)的基數(shù)為b的冪增型泛水仙花數(shù),且顯然位數(shù)為1時(shí)只有1是泛水仙花數(shù)，故不再計(jì)算位數(shù)為1的情況。

給定位數(shù)n，計(jì)算8位數(shù)以內(nèi)的所有冪增型泛水仙花數(shù)需要7141.182828秒即約為119分鐘的時(shí)間。表3列出了108以內(nèi)的冪等差水仙花數(shù)。

表3 108以內(nèi)的冪等差數(shù)

3 結(jié)語

本文將林宣治、楊健等研究的齊次水仙花數(shù)推廣為更一般的冪等差型水仙花數(shù)，使得齊次水仙花數(shù)成為我們研究范疇的一個(gè)特例，并對冪等差型水仙花數(shù)的性質(zhì)進(jìn)行了研究。

在對冪等差型水仙花數(shù)研究的過程中，不禁產(chǎn)生了一些疑惑，其中最大的疑惑就是：

問題1是否存在任意公差的冪等差型水仙花數(shù)？

顯然，無論計(jì)算機(jī)技術(shù)發(fā)展到什么程度都無法從計(jì)算的角度解決上述問題。于是，我們把上述問題等價(jià)地轉(zhuǎn)化為如下數(shù)學(xué)問題：

問題2下述關(guān)于(n;b,d,x1,x2,…,xn)的不定方程

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

期待數(shù)學(xué)工作者給出問題2的解答。