水仙花数

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}}}