【讨论】关于水仙花程序的matlab设计

qibbxxt

17# shamohu
恩，lingo和1stOpt应该可以做，但是matlab做不了整数规划，而去设计演化算法对这一类问题来说，比较难，因为很难定义解的领域，来选择或者构造相应的操作算子，不过启发式算法一旦设计好了，对于大规模的这一类问题应该就比较方便了

白水湾湾

插一句嘴，不知道对不对，楼主为什么喜欢用fliplr这个函数呢？？
clear;clc;close all
tic
a=100:999;b=0:2;
a(a'==sum([mod(floor(bsxfun(@rdivide,[a' a' a'],10.^fliplr(b))),10)].^3,2))
toc

既然是追求效率，直接把b赋值成2:-1:0效率不更高吗？？

qibbxxt

22# 白水湾湾
谢谢，起初这是为了得到结果
后来就忘记了这个问题

qibbxxt

17# shamohu
我用matlab的最优化方法去求解，但是算法容易陷入局部最小

1. 
   
2. function y=myfun(x)
   
3. n=3;
   
4. b=n-1:-1:0;
   
5. y=(sum(mod(floor(bsxfun(@rdivide,repmat(x,1,n),10.^fliplr(b))),10).^n,2)-x).^2;
   

复制代码

用下面的

1. 
   
2. >> [x,y]=fminsearch(@myfun,400)
   
3. x =
   
4.    370
   
5. 
   
6. y =
   
7.      0
   

复制代码

针对多个局部极值，我用遗传工具箱，但是位数多的时候就不好用了

1. 
   
2. function [x,fval,exitflag,output,population,score] = myfile(nvars,lb,ub,PopulationSize_Data)
   
3. % This is an auto generated M-file from Optimization Tool.
   
4. % Start with the default options
   
5. options = gaoptimset;
   
6. % Modify options setting
   
7. options = gaoptimset(options,'PopulationSize', PopulationSize_Data);
   
8. options = gaoptimset(options,'Display', 'off');
   
9. [x,fval,exitflag,output,population,score] = ...
   
10. ga(@myfun,nvars,[],[],[],[],lb,ub,[],options);
    

复制代码

lin2009

“乐在其中的数学”提到国防科技大学的刘江宁用了18小时算出全部的水仙花数（书中的页码282～286）。

Lz看看。

注： “乐在其中的数学”文件太大无法上传，请用百度搜索下载。

qibbxxt

25# lin2009
谢谢，有机会我下载下来看一下，顺便看看他是基于什么思想编写的程序

qibbxxt

25# lin2009
由于计算机内存问题，向量化内存就会报错，而且repmat的效率太低，我将16#的程序该作循环，空间化为时间，运算了一下，7位数如下：

1. 
   
2. clear;clc
   
3. n=7;
   
4. tic
   
5. a=10^(n-1):10^n-1;
   
6. b=(0:9).^n;
   
7. for i=1:length(a)
   
8.     if sum(b(mod(floor(a(i)./10.^((n-1):-1:0)),10)+1))==a(i)
   
9.         disp(a(i));
   
10.     end
    
11. end
    
12. toc
    

复制代码

1. 
   
2.      1741725
   
3.      4210818
   
4.      9800817
   
5.      9926315
   
6. Elapsed time is 41.598189 seconds.
   

复制代码

lin2009

27# qibbxxt
分解为7 = 1 + 6的问题可能更快，前面的1，即第7位用循环，后面6位用矩阵／数组操作。

你的机子比我快。作为对比应该是在同一台电脑上运行同一程序。
你的程序在我的机子上运行了295秒。（你的41秒，呵呵）
     1741725
     4210818
     9800817
     9926315
Elapsed time is 295.893893 seconds.

下面的程序（7 = 1 + 6程序）在我的机子上运行的时间为36.275952（全部都找到了，但是后面输出的格式有些不对，有点莫名其妙）。
在你的机子上应该更快！

1. 
   
2. tic
   
3. [d1,d2,d3,d4,d5,d6]  = ndgrid(0:9,0:9,0:9,0:9,0:9,0:9);
   
4. for d0 = 1:9
   
5.     ind = ( d0^7 + d1.^7 + d2.^7 + d3.^7 + d4.^7 + d5.^7 + d6.^7 == 10^6*d0 + 10^5*d1 + 10^4*d2 + 10^3*d3 + 10^2*d4 +10^1*d5 + d6);   
   
6.     Num = 10^6*d0 + 10^5*d1(ind) + 10^4*d2(ind) + 10^3*d3(ind) + 10^2*d4(ind) + 10^1*d5(ind) +d6(ind);
   
7.     if ~isempty(Num)
   
8.         fprintf('\t%d = %d^7 + %d^7 + %d^7 + %d^7 + %d^7 + %d^7 + %d^7\n',Num,d0,d1(ind),d2(ind),d3(ind),d4(ind),d5(ind),d6(ind));   
   
9.     end
   
10. end
    
11. toc
    

复制代码

1741725 = 1^7 + 7^7 + 4^7 + 1^7 + 7^7 + 2^7 + 5^7
4210818 = 4^7 + 2^7 + 1^7 + 0^7 + 8^7 + 1^7 + 8^7
9926315 = 9800817^7 + 9^7 + 9^7 + 8^7 + 2^7 + 0^7 + 6^7
0 = 3^7 + 8^7 + 1^7 + 1^7 + 5^7 + 7^7 + Elapsed time is 36.275952 seconds.

qibbxxt

28# lin2009
的确这样计算的速度是要快一些，你的程序在我的电脑上面运行结果如下：

1. 
   
2. 1741725 = 1^7 + 7^7 + 4^7 + 1^7 + 7^7 + 2^7 + 5^7
   
3. 4210818 = 4^7 + 2^7 + 1^7 + 0^7 + 8^7 + 1^7 + 8^7
   
4. 9926315 = 9800817^7 + 9^7 + 9^7 + 8^7 + 2^7 + 0^7 + 6^7
   
5. 0 = 3^7 + 8^7 + 1^7 + 1^7 + 5^7 + 7^7 + Elapsed time is 8.271417 seconds.
   

复制代码

qibbxxt

28# lin2009
如果直接按照我的第二种思路写的话，要慢一些，看来分解会加快速度

1. 
   
2. clear;clc;close all
   
3. tic;
   
4. [x1,x2,x3,x4,x5,x6,x7]=ndgrid(1:9,0:9,0:9,0:9,0:9,0:9,0:9);
   
5. ind=(x1.^7+x2.^7+x3.^7+x4.^7++x5.^7+x6.^7+x7.^7==1000000*x1+100000*x2+10000*x3+1000*x4+100*x5+10*x6+x7);
   
6. 1000000*x1(ind)+100000*x2(ind)+10000*x3(ind)+1000*x4(ind)+100*x5(ind)+10*x6(ind)+x7(ind)
   
7. toc
   

复制代码

1. 
   
2. ans =
   
3.      9926315
   
4.      1741725
   
5.      9800817
   
6.      4210818
   
7. Elapsed time is 10.822329 seconds.
   

复制代码

qibbxxt

28# lin2009
你的程序问题在于以9开头的水仙花数有2个，所以就乱了
我把leng版的思想和我的第二种思想结合，写出以下2个程序，分别是

1. 
   
2. clear;clc;close all
   
3. tic;
   
4. b=(0:9).^7;
   
5. [x1,x2,x3,x4,x5,x6,x7]=ndgrid(1:9,0:9,0:9,0:9,0:9,0:9,0:9);
   
6. x=b(cat(8,x1,x2,x3,x4,x5,x6,x7)+1);
   
7. y=sum(x,8);
   
8. ind=(y==1000000*x1+100000*x2+10000*x3+1000*x4+100*x5+10*x6+x7);
   
9. 1000000*x1(ind)+100000*x2(ind)+10000*x3(ind)+1000*x4(ind)+100*x5(ind)+10*x6(ind)+x7(ind)
   
10. toc
    

复制代码

1. 
   
2. ans =
   
3.      9926315
   
4.      1741725
   
5.      9800817
   
6.      4210818
   
7. Elapsed time is 4.519278 seconds.
   

复制代码

和

1. 
   
2. clear;clc;close all
   
3. tic;
   
4. b=(0:9).^7;
   
5. [x1,x2,x3,x4,x5,x6,x7]=ndgrid(1:9,0:9,0:9,0:9,0:9,0:9,0:9);
   
6. ind=(b(x1+1)+b(x2+1)+b(x3+1)+b(x4+1)++b(x5+1)+b(x6+1)+b(x7+1)==1000000*x1+100000*x2+10000*x3+1000*x4+100*x5+10*x6+x7);
   
7. 1000000*x1(ind)+100000*x2(ind)+10000*x3(ind)+1000*x4(ind)+100*x5(ind)+10*x6(ind)+x7(ind)
   
8. toc
   

复制代码

1. 
   
2. 
   
3. ans =
   
4.      9926315
   
5.      1741725
   
6.      9800817
   
7.      4210818
   
8. Elapsed time is 4.331290 seconds.
   

复制代码

shamohu

在下面链接发了个类似的趣味题，有兴趣的看看：
http://forum.simwe.com/viewthread.php?tid=949103&extra=