关于完全数的matlab程序设计

Everley

来自：百度百科简介　　完全数（Perfect number），又称完美数或完备数，是一些特殊的自然数：
它所有的真因子（即除了自身以外的约数）的和（即因子函数），恰好等于它本身。
举例　　例如：第一个完全数是6，它有约数1、2、3、6，除去它本身6外，其余3个数相加，1+2+3＝6。第二个完全数是28，它有约数1、2、4、7、14、28，除去它本身28外，其余5个数相加，1+2+4+7+14＝28。后面的数是496、8128等等。

　　例如：

　　6＝1+2+3

　　28＝1+2+4+7+14

　　496＝1+2+4+8+16+31+62+124+248

　　8128＝1+2+4+8+16+32+64+127+254+508+1016+2032+4064
那如何用matlab找到这些数字呢？
1.能否少做循环或者不用循环
2.如何提高程序运行效率
以下是我的做法，希望大家一起讨论
第一种：

1. clear;clc;close all
   
2. a=1:10000;
   
3. a(2*a'==accumarray(a',a',[],@(x)sum(a(abs(x./(1:x)-round(x./(1:x)))<eps))));
   

复制代码

第二种：

1. a=1:14;
   
2. c=(2.^a(isprime(2.^a-1))-1).*2.^(a(isprime(2.^a-1))-1)
   

复制代码

第三种：

1. 
   
2. a=1:10000;
   
3. d=a(2*a'==accumarray(a',a',[],@(x)prod((bsxfun(@power,find(accumarray(factor(x)',1)),feval(@(y)y(y>0),accumarray(factor(x)',1))+1)-1)./(find(accumarray(factor(x)',1))-1))));
   

复制代码

第四种：

1. a(2*a'==accumarray(a',a',[],@(x)sum(a(mod(x,1:x)==0))))

复制代码

Claudia

1. 
   
2. tic
   
3. n=10000;
   
4. for i=1:n
   
5. if sum(find(bsxfun(@mod,i,1:i-1)==0))==i
   
6. disp(i)
   
7. end
   
8. end
   
9. toc
   

复制代码

6

    28

   496

        8128

Elapsed time is 7.991244 seconds.       我的电脑太慢了。。 .循环液尽量少用，费时！~

Amaris

解法二，用的是完整数的公式。
受他启发，直接google查询最快，呵呵。

Euclid discovered that the first four perfect numbers are generated by the formula

                               
登录/注册后可看大图

:
....
The first 39 even perfect numbers are 2^(n&#8722;1) * (2^n &#8722; 1) for
n = 2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917
n = [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917];
perfectNo = 2.^(n-1).*(2.^n - 1);

第 13 个完整数(n=521，对应的完整数为2^521*(2^250-1)，就超出计算的运算范围！（Matlab显示 Inf）

Everley

上次在想如何去掉数组中的0元素，又不影响数组作为整体参数的操作，所以用

1. feval(@(y)y(y>0),accumarray(factor(x)',1))+1)

复制代码

今天看到了nonzeros,所以上面的语句就可以替代为

1. nonzeros(accumarray(factor(x)',1))+1)

复制代码

整体就会缩短一点，为

1. clear
   
2. a=1:10000;
   
3. d=a(2*a'==accumarray(a',a',[],@(x)prod((bsxfun(@power,find(accumarray(factor(x)',1)),nonzeros(accumarray(factor(x)',1))+1)-1)./(find(accumarray(factor(x)',1))-1))));

复制代码

Amaris

与 完整数 相关的是 亲和数。

用计算机找出10^6(即100万以内的亲和数对)，看看时间花多少？

1. 
   
2. clear all;
   
3. clc;
   
4. start_time = tic;
   
5. fid1 = fopen('wanquanshu.m', 'w'); % 完全数
   
6. fid2 = fopen('qinheshu.m', 'w'); % 亲和数
   
7. for n = 10:10^6
   
8.     if isprime(n)
   
9.         continue;
   
10.     end
    
11.     m = sum(find(fix(n ./ (1 : n)) == n ./ (1 : n))) - n; % 约数之和, 低效的算法！
    
12.     if m == n
    
13.         fprintf(fid1, ' %d\n', n);
    
14.         continue;
    
15.     end
    
16.     nn = sum(find(fix(m ./ (1 : m)) == m ./ (1 : m))) - m;
    
17.     if nn == n
    
18.         fprintf(fid2, ' %d, %d\n', n, m);
    
19.     end
    
20. end
    
21. time_comsuming = toc;
    
22. fprintf(fid1, '\n Elapsed time is %f seconds.\n', time_comsuming);
    
23. fprintf(fid1, '\n Elapsed time is %f seconds.\n', time_comsuming);
    
24. fclose(fid1);
    
25. fclose(fid2);
    
26. !shutdown - s
    

复制代码

在Celeron(R) CPU 2.66G 760MB内存上，耗时超过24小时！（有点无聊，呵呵）

1、大家看看有没有高效的算法？应该可以利用factor(n)及约数的通式a^n1*b^n2*c^n3...以及accumarray及bsxfun等函数？

结果如下，
2、怎样把重复的亲和数对去掉？

1. 220, 284
   
2. 284, 220
   
3. 1184, 1210
   
4. 1210, 1184
   
5. 2620, 2924
   
6. 2924, 2620
   
7. 5020, 5564
   
8. 5564, 5020
   
9. 6232, 6368
   
10. 6368, 6232
    
11. 10744, 10856
    
12. 10856, 10744
    
13. 12285, 14595
    
14. 14595, 12285
    
15. 17296, 18416
    
16. 18416, 17296
    
17. 63020, 76084
    
18. 66928, 66992
    
19. 66992, 66928
    
20. 67095, 71145
    
21. 69615, 87633
    
22. 71145, 67095
    
23. 76084, 63020
    
24. 79750, 88730
    
25. 87633, 69615
    
26. 88730, 79750
    
27. 100485, 124155
    
28. 122265, 139815
    
29. 122368, 123152
    
30. 123152, 122368
    
31. 124155, 100485
    
32. 139815, 122265
    
33. 141664, 153176
    
34. 142310, 168730
    
35. 153176, 141664
    
36. 168730, 142310
    
37. 171856, 176336
    
38. 176272, 180848
    
39. 176336, 171856
    
40. 180848, 176272
    
41. 185368, 203432
    
42. 196724, 202444
    
43. 202444, 196724
    
44. 203432, 185368
    
45. 280540, 365084
    
46. 308620, 389924
    
47. 319550, 430402
    
48. 356408, 399592
    
49. 365084, 280540
    
50. 389924, 308620
    
51. 399592, 356408
    
52. 430402, 319550
    
53. 437456, 455344
    
54. 455344, 437456
    
55. 469028, 486178
    
56. 486178, 469028
    
57. 503056, 514736
    
58. 514736, 503056
    
59. 522405, 525915
    
60. 525915, 522405
    
61. 600392, 669688
    
62. 609928, 686072
    
63. 624184, 691256
    
64. 635624, 712216
    
65. 643336, 652664
    
66. 652664, 643336
    
67. 667964, 783556
    
68. 669688, 600392
    
69. 686072, 609928
    
70. 691256, 624184
    
71. 712216, 635624
    
72. 726104, 796696
    
73. 783556, 667964
    
74. 796696, 726104
    
75. 802725, 863835
    
76. 863835, 802725
    
77. 879712, 901424
    
78. 898216, 980984
    
79. 901424, 879712
    
80. 947835, 1125765
    
81. 980984, 898216
    
82. 998104, 1043096

复制代码

Everley

1.请看我关于完全数的第三种做法，就是用向量操作计算约数和的另一种方法，就是用factor实现的，基于完全数的程序，加几行代码就可以实现亲和数
2.只需要取所得数的奇数行或者偶数行
因此我的代码如下

1. 
   
2. clear;clc;close all
   
3. clear
   
4. a=1:10000;
   
5. d=accumarray(a',a',[],@(x)prod((bsxfun(@power,find(accumarray(factor(x)',1)),nonzeros(accumarray(factor(x)',1))+1)-1)./(find(accumarray(factor(x)',1))-1)))-a';
   
6. e=ones(max(d),1);
   
7. e(a)=d;e(1)=1;
   
8. f=[e(e(e(a))==a'),e(e(e(e(a))==a'))];
   
9. f(f(:,1)==f(:,2),:)=[];
   
10. g=f(2:2:end,:)
    

复制代码

1. 
   
2. g =
   
3.          220         284
   
4.         1184        1210
   
5.         2620        2924
   
6.         5020        5564
   
7.         6232        6368
   
8. 
   

复制代码

Everley

1. g =
   
2. 
   
3.          220         284
   
4.         1184        1210
   
5.         2620        2924
   
6.         5020        5564
   
7.         6232        6368
   
8.        10744       10856
   
9.        12285       14595
   
10.        17296       18416
    
11.        66992       66928
    
12.        71145       67095
    
13.        67095       71145
    
14.        88730       79750
    
15.        79750       88730
    
16.       139815      122265
    
17.       122368      123152
    
18.       122265      139815
    
19.       168730      142310
    
20.       142310      168730
    
21.       180848      176272
    
22.       176272      180848
    
23.       202444      196724
    
24.       185368      203432
    
25.       389924      308620
    
26.       399592      356408
    
27.       308620      389924
    
28.       319550      430402
    
29.       437456      455344
    
30.       469028      486178
    
31.       503056      514736
    
32.       522405      525915
    
33.       686072      609928
    
34.       712216      635624
    
35.       643336      652664
    
36.       600392      669688
    
37.       624184      691256
    
38.       796696      726104
    
39.       726104      796696
    
40.       802725      863835
    
41.       980984      898216
    
42.       898216      980984
    
43. 
    
44. Elapsed time is 1217.180373 seconds.

复制代码

这是在我的电脑上运行1000000内的结果

Amaris

Everley 发表于 2016-9-14 11:25
1.请看我关于完全数的第三种做法，就是用向量操作计算约数和的另一种方法，就是用factor实现的，基于完全数 ...

代码简洁，效率高，赞同rocwoods的说法：“齐老弟的代码让我想起 GRE考试、各地方言水平能力测试题来了。如果MATLAB也有变态的水平测试题的话，上面代码可以入选啦”。


不过第二问，不是简单的“只需要取所得数的奇数行或者偶数行”
前面几行是这样的，但是后面的没有这样的规律。

如7#的25～28行,等

87633, 69615
88730, 79750
100485, 124155
122265, 139815

Everley

纠正一下，取奇数或者偶数行没有办法解决问题2
所以代码改为

1. 
   
2. clear;clc;close all
   
3. clear
   
4. a=1:10000;
   
5. d=accumarray(a',a',[],@(x)prod((bsxfun(@power,find(accumarray(factor(x)',1)),nonzeros(accumarray(factor(x)',1))+1)-1)./(find(accumarray(factor(x)',1))-1)))-a';
   
6. e=ones(max(d),1);
   
7. e(a)=d;e(1)=1;
   
8. f=[e(e(e(a))==a'),e(e(e(e(a))==a'))];
   
9. f(f(:,1)==f(:,2),:)=[];
   
10. h=sum(f,2);
    
11. [n, m] = unique(h, 'first');
    
12. p=f(sort(m),:)
    

复制代码

得

1. p =
   
2. 
   
3.          284         220
   
4.         1210        1184
   
5.         2924        2620
   
6.         5564        5020
   
7.         6368        6232
   
8.        10856       10744
   
9.        14595       12285
   
10.        18416       17296
    
11.        76084       63020
    
12.        66992       66928
    
13.        71145       67095
    
14.        87633       69615
    
15.        88730       79750
    
16.       124155      100485
    
17.       139815      122265
    
18.       123152      122368
    
19.       153176      141664
    
20.       168730      142310
    
21.       176336      171856
    
22.       180848      176272
    
23.       203432      185368
    
24.       202444      196724
    
25.       365084      280540
    
26.       389924      308620
    
27.       430402      319550
    
28.       399592      356408
    
29.       455344      437456
    
30.       486178      469028
    
31.       514736      503056
    
32.       525915      522405
    
33.       669688      600392
    
34.       686072      609928
    
35.       691256      624184
    
36.       712216      635624
    
37.       783556      667964
    
38.       796696      726104
    
39.       863835      802725
    
40.       901424      879712
    
41.       980984      898216

复制代码

Everley

一种更好的处理办法如下

1. 
   
2. >> q=unique(sort(f,2),'rows')
   
3. q =
   
4.          220         284
   
5.         1184        1210
   
6.         2620        2924
   
7.         5020        5564
   
8.         6232        6368
   
9.        10744       10856
   
10.        12285       14595
    
11.        17296       18416
    
12.        63020       76084
    
13.        66928       66992
    
14.        67095       71145
    
15.        69615       87633
    
16.        79750       88730
    
17.       100485      124155
    
18.       122265      139815
    
19.       122368      123152
    
20.       141664      153176
    
21.       142310      168730
    
22.       171856      176336
    
23.       176272      180848
    
24.       185368      203432
    
25.       196724      202444
    
26.       280540      365084
    
27.       308620      389924
    
28.       319550      430402
    
29.       356408      399592
    
30.       437456      455344
    
31.       469028      486178
    
32.       503056      514736
    
33.       522405      525915
    
34.       600392      669688
    
35.       609928      686072
    
36.       624184      691256
    
37.       635624      712216
    
38.       643336      652664
    
39.       667964      783556
    
40.       726104      796696
    
41.       802725      863835
    
42.       879712      901424
    
43.       898216      980984
    

复制代码

代码中的f就是7#的f

陌影

Amaris 发表于 2016-9-14 11:25
与 完整数 相关的是 亲和数。

用计算机找出10^6(即100万以内的亲和数对)，看看时间花多少？

耗时24小时  哈哈哈  费计算机啊

陌影

Everley 发表于 2016-9-14 11:28
纠正一下，取奇数或者偶数行没有办法解决问题2
所以代码改为
得

学习了  好方法