征求“水仙花数”高效算法

chyanog

可能大家都知道，绝大多数编程基础书中都有水仙花数这一经典问题 。广义的水仙花数是指一个 n 位数 ( n≥3 )，它的每个位上的数字的 n 次幂之和等于它本身。（例如：1^3 + 5^3 + 3^3 = 153），它的求法可真不少。我用M@写了几种，先以三位的为例，如下：
（1）

1. 
   
2. For[a = 1, a <= 9, a++,
   
3. For[b = 0, b <= 9, b++,
   
4.   For[c = 0, c <= 9, c++,
   
5.    If[FromDigits[{a, b, c}] == Total[{a, b, c}^3], Print[a, b, c]]]]]

复制代码

（2）

1. 
   
2. For[n = 100, n < 1000, n++,
   
3.   a = Mod[n, 10];(*IntegerDigits[n][[1]]*)
   
4.   b = Mod[IntegerPart[n/10], 10]; (*Mod[Quotient[i,10], 10] 或IntegerDigits[n][[2]]*)
   
5.   c = IntegerPart[n/100];(*Quotient[i,100]或IntegerDigits[n][[3]]*)
   
6.   If[a^3 + b^3 + c^3 == n,
   
7.    Print[n]]] // Timing

复制代码

（3）

1. 
   
2. Select[Tuples[Range[0, 9], 3],
   
3.   Tr[{#[[1]], #[[2]], #[[3]]}^3] ==
   
4.     FromDigits[{#[[1]], #[[2]], #[[3]]}] &] // Timing

复制代码

（4）

1. 
   
2. Reduce[a^3 + b^3 + c^3 == FromDigits[{a, b, c}] && 1 <= a <= 9 &&
   
3.    0 <= b <= 9 && 0 <= c <= 9, {a, b, c}, Integers] // Timing

复制代码

以下两个为不限位数的算法：
（5）

1. 
   
2. Module[{n, i, k, s, t1, t2},
   
3.   n = Input["输入不小于3的自然数n:"];
   
4.   i = 10^(n - 1);
   
5.   While[i <= 10^n - 1,
   
6.    s = 0;
   
7.    k = 1;
   
8.    While[k <= n,
   
9.     t1 = IntegerPart[i/10^(k - 1)]; (*Floor[i/10^(k-1)];*)
   
10.     t2 = Mod[t1, 10]^n;
    
11.     s = s + t2;
    
12.     k++];
    
13.    If[s == i, Print];
    
14. i++];
    
15.   ] // Timing

复制代码

(6)

1. 
   
2.   Module[{n},
   
3.   n = Input["请输入不小于3的自然数:"];
   
4.   Select[Range[10^(n - 1), 10^n - 1],
   
5.    Tr[IntegerDigits[#]^n] == # &]] // Timing

复制代码

其中最后一个是我写的效率相对较高也最简洁的啦，（之所以用Tr是因为多次测试发现它往往比Total更快，有点儿郁闷）
不过，当所求的水仙花数位数为7时，运行时间就几乎60s了，位数更高时就难以想象了。不知道我的代码是否可以再优化一下，或者还有其他精妙的算法，恳请诸位方家赐教！

ggggwhw

并不是每个人都知道所谓的水仙花数的概念的,
尤其是你所谓的位数不等于3时的水仙花数的概念.
为什么你的第6种情况下要求输入的数字不小于3呢?

chyanog

3# waynebuaa
不知这位朋友用的是哪种语言，算法上可否告知呢？

liweicai990

1. n = Input["Input  >2  integer n:"];
   
2. For[i = 10^(n - 1), i < 10^n - 1, i++,
   
3. If[Total[IntegerDigits[i]^IntegerLength[i]] == i, Print[i]]]

复制代码