水仙花数

sir_chen

水仙花数原定义是这样的:对于一个三位数,若各位数的立方和等于这个数本身,那么这个数就是水仙花数.
现在对这个定义做一个扩展,即对于一个n位数,若各位数的n次方和等于这个数本身,那么这个n位数就是一个水仙花数.
例如:对于4位数$1643=1^4+6^4+3^4+4^4$,从而1634是一个4位水仙花数.
现在有3个问题:
(1)最大存在多少位的水仙花数
(2)最大的水仙花数是多少
(3)一共有多少个水仙花数

sir_chen

对于问题1,最大存在59位的水仙花数
用$a_k$表示每个数位,若一个数是水仙花数则满足:
$sum_{k=0}^{n}a_k r^k=sum_{k=0}^{n}a_k^(n+1)$(此处r=10是基数)
也即$sum_{k=0}^{n}(a_k^(n+1)-a_k r^k)=0$
记$a=(a_0,a_1,......,a_n)$
$f(n,a)=sum_{k=0}^{n}(a_k^(n+1)-a_k r^k)$
对于一个n+1位数,$a_n>=1$
当$a_k=0,(k<=n-1)$时,有$f(n,a)=a_n^(n+1)-a_n r^n=a_n(a_n^n-r^n)<0$
因此只有当f(n,a)的最大值不小于0时,f(n,a)=0才可能有根
考察函数$g(x)=x^(n+1)-x r^k,0<=x<=r-1$,此函数在$x=r^(k/n)/root{n}{n+1}$处取到最小值,最大值在端点处取得.
所以$g_max(x) =y(k)={(0,k>{n ln(r-1)}/lnr),((r-1)((r-1)^n-r^k),k<={n ln(r-1)}/lnr):}$
当k=n时,由于$a_n>=1$,所以$y(n)=1-r^n$
所以$f_max(n,a)=1-r^n+sum_{k=0}^(n-1)y(k)$
当$f_max(n,a)<0 $时,则方程f(n,a)=0没有根
编程求解得到$n<=59$
所以最大存在60位的水仙花数

sir_chen

对于1-9位的水仙花数如下
1
2
3
4
5
6
7
8
9
153
370
371
407
1634
8208
9474
54748
92727
93084
548834
1741725
4210818
9800817
9926315
24678050
24678051
88593477
146511208
472335975
534494836
912985153
由于程序采用的是穷举法,所以到10位以后运行较慢,暂无满意结果

wiley

A005188

10进制下共有88个:
http://www.research.att.com/~njas/sequences/b005188.txt

or:
Narcissistic Number

sir_chen

只是简单的列出了水仙花数,但是没有说是怎么找出来的

zed

http://bbs.emath.ac.cn/viewthrea ... ;fromuid=41#pid1427

mathematica

Clear["Global`*"];(*Clear all variables*)
Select[Range[1,5000],#==Total@(IntegerDigits[#]^3)&]
这个是代码

mathematica

Clear["Global`*"];(*Clear all variables*)
Select[Range[1,500000],#==Total@(IntegerDigits[#]^(Length@IntegerDigits[#]))&]
{1, 2, 3, 4, 5, 6, 7, 8, 9, 153, 370, 371, 407, 1634, 8208, 9474, \
54748, 92727, 93084} 这个是运算结果

mathematica

英语翻译是happy number

chyanog

mathematica 发表于 2013-10-14 13:05
Clear["Global`*"];(*Clear all variables*)
Select[Range[1,500000],#==Total@(IntegerDigits[#]^(Length ...

这样比较慢，可以用Compile加速（不用Compile的代码在我的电脑上5s）

1. Compile[{},
   
2.    Select[Range[1, 500000],
   
3.     With[{id = IntegerDigits@#}, # == Total[id^Length@id]] &]
   
4.    ][] // AbsoluteTiming

复制代码

{0.288017, {1, 2, 3, 4, 5, 6, 7, 8, 9, 153, 370, 371, 407, 1634, 8208, 9474, 54748, 92727, 93084}}

试试这个（7位以内的只需0.3s）：

1. nar[n_] := Block[{A},
   
2.    A = Flatten@Outer[Plus, ##] & @@ Array[Range[Boole[# == 1], 9]^n &, n];
   
3.    Pick[A, A~BitXor~Range[10^(n - 1), 10^n - 1], 0]
   
4.    ];
   
5. nar /@ Range[3, 7] // AbsoluteTiming

复制代码

{0.336019, {{153, 370, 371, 407}, {1634, 8208, 9474}, {54748, 92727,  93084}, {548834}, {1741725, 4210818, 9800817, 9926315}}}