关于完全数的matlab程序设计
来自：百度百科简介 完全数（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.如何提高程序运行效率
以下是我的做法，希望大家一起讨论
第一种：
clear;clc;close all
a=1:10000;
a(2*a'==accumarray(a',a',[],@(x)sum(a(abs(x./(1:x)-round(x./(1:x)))<eps))));
第二种：
a=1:14;
c=(2.^a(isprime(2.^a-1))-1).*2.^(a(isprime(2.^a-1))-1)
第三种：
a=1:10000;
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))));
第四种：
a(2*a'==accumarray(a',a',[],@(x)sum(a(mod(x,1:x)==0))))
tic
n=10000;
for i=1:n
if sum(find(bsxfun(@mod,i,1:i-1)==0))==i
disp(i)
end
end
toc
6
28
496
8128
Elapsed time is 7.991244 seconds. 我的电脑太慢了。。 .循环液尽量少用，费时！~
解法二，用的是完整数的公式。
受他启发，直接google查询最快，呵呵。
Euclid discovered that the first four perfect numbers are generated by the formulahttp://latex.codecogs.com/gif.latex?${{2}^{n-1}}\left(%20{{2}^{n}}-1%20\right)$:
....
The first 39 even perfect numbers are 2^(n−1) * (2^n − 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 = ;
perfectNo = 2.^(n-1).*(2.^n - 1);
第 13 个完整数(n=521，对应的完整数为2^521*(2^250-1)，就超出计算的运算范围！（Matlab显示 Inf）
上次在想如何去掉数组中的0元素，又不影响数组作为整体参数的操作，所以用
feval(@(y)y(y>0),accumarray(factor(x)',1))+1)今天看到了nonzeros,所以上面的语句就可以替代为
nonzeros(accumarray(factor(x)',1))+1)整体就会缩短一点，为
clear
a=1:10000;
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))));
与 完整数 相关的是 亲和数。
用计算机找出10^6(即100万以内的亲和数对)，看看时间花多少？
clear all;
clc;
start_time = tic;
fid1 = fopen('wanquanshu.m', 'w'); % 完全数
fid2 = fopen('qinheshu.m', 'w'); % 亲和数
for n = 10:10^6
if isprime(n)
continue;
end
m = sum(find(fix(n ./ (1 : n)) == n ./ (1 : n))) - n; % 约数之和, 低效的算法！
if m == n
fprintf(fid1, ' %d\n', n);
continue;
end
nn = sum(find(fix(m ./ (1 : m)) == m ./ (1 : m))) - m;
if nn == n
fprintf(fid2, ' %d, %d\n', n, m);
end
end
time_comsuming = toc;
fprintf(fid1, '\n Elapsed time is %f seconds.\n', time_comsuming);
fprintf(fid1, '\n Elapsed time is %f seconds.\n', time_comsuming);
fclose(fid1);
fclose(fid2);
!shutdown - s
在Celeron(R) CPU 2.66G 760MB内存上，耗时超过24小时！（有点无聊，呵呵）
1、大家看看有没有高效的算法？应该可以利用factor(n)及约数的通式a^n1*b^n2*c^n3...以及accumarray及bsxfun等函数？
结果如下，
2、怎样把重复的亲和数对去掉？
220, 284
284, 220
1184, 1210
1210, 1184
2620, 2924
2924, 2620
5020, 5564
5564, 5020
6232, 6368
6368, 6232
10744, 10856
10856, 10744
12285, 14595
14595, 12285
17296, 18416
18416, 17296
63020, 76084
66928, 66992
66992, 66928
67095, 71145
69615, 87633
71145, 67095
76084, 63020
79750, 88730
87633, 69615
88730, 79750
100485, 124155
122265, 139815
122368, 123152
123152, 122368
124155, 100485
139815, 122265
141664, 153176
142310, 168730
153176, 141664
168730, 142310
171856, 176336
176272, 180848
176336, 171856
180848, 176272
185368, 203432
196724, 202444
202444, 196724
203432, 185368
280540, 365084
308620, 389924
319550, 430402
356408, 399592
365084, 280540
389924, 308620
399592, 356408
430402, 319550
437456, 455344
455344, 437456
469028, 486178
486178, 469028
503056, 514736
514736, 503056
522405, 525915
525915, 522405
600392, 669688
609928, 686072
624184, 691256
635624, 712216
643336, 652664
652664, 643336
667964, 783556
669688, 600392
686072, 609928
691256, 624184
712216, 635624
726104, 796696
783556, 667964
796696, 726104
802725, 863835
863835, 802725
879712, 901424
898216, 980984
901424, 879712
947835, 1125765
980984, 898216
998104, 1043096
1.请看我关于完全数的第三种做法，就是用向量操作计算约数和的另一种方法，就是用factor实现的，基于完全数的程序，加几行代码就可以实现亲和数
2.只需要取所得数的奇数行或者偶数行
因此我的代码如下
clear;clc;close all
clear
a=1:10000;
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';
e=ones(max(d),1);
e(a)=d;e(1)=1;
f=;
f(f(:,1)==f(:,2),:)=[];
g=f(2:2:end,:)
g =
220 284
1184 1210
2620 2924
5020 5564
6232 6368
g =
220 284
1184 1210
2620 2924
5020 5564
6232 6368
10744 10856
12285 14595
17296 18416
66992 66928
71145 67095
67095 71145
88730 79750
79750 88730
139815 122265
122368 123152
122265 139815
168730 142310
142310 168730
180848 176272
176272 180848
202444 196724
185368 203432
389924 308620
399592 356408
308620 389924
319550 430402
437456 455344
469028 486178
503056 514736
522405 525915
686072 609928
712216 635624
643336 652664
600392 669688
624184 691256
796696 726104
726104 796696
802725 863835
980984 898216
898216 980984
Elapsed time is 1217.180373 seconds.这是在我的电脑上运行1000000内的结果
Everley 发表于 2016-9-14 11:25
1.请看我关于完全数的第三种做法，就是用向量操作计算约数和的另一种方法，就是用factor实现的，基于完全数 ...
代码简洁，效率高，赞同rocwoods的说法：“齐老弟的代码让我想起 GRE考试、各地方言水平能力测试题来了。如果MATLAB也有变态的水平测试题的话，上面代码可以入选啦”。
不过第二问，不是简单的“只需要取所得数的奇数行或者偶数行”
前面几行是这样的，但是后面的没有这样的规律。
如7#的25～28行,等
87633, 69615
88730, 79750
100485, 124155
122265, 139815 纠正一下，取奇数或者偶数行没有办法解决问题2
所以代码改为
clear;clc;close all
clear
a=1:10000;
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';
e=ones(max(d),1);
e(a)=d;e(1)=1;
f=;
f(f(:,1)==f(:,2),:)=[];
h=sum(f,2);
= unique(h, 'first');
p=f(sort(m),:)
得
p =
284 220
1210 1184
2924 2620
5564 5020
6368 6232
10856 10744
14595 12285
18416 17296
76084 63020
66992 66928
71145 67095
87633 69615
88730 79750
124155 100485
139815 122265
123152 122368
153176 141664
168730 142310
176336 171856
180848 176272
203432 185368
202444 196724
365084 280540
389924 308620
430402 319550
399592 356408
455344 437456
486178 469028
514736 503056
525915 522405
669688 600392
686072 609928
691256 624184
712216 635624
783556 667964
796696 726104
863835 802725
901424 879712
980984 898216
一种更好的处理办法如下
>> q=unique(sort(f,2),'rows')
q =
220 284
1184 1210
2620 2924
5020 5564
6232 6368
10744 10856
12285 14595
17296 18416
63020 76084
66928 66992
67095 71145
69615 87633
79750 88730
100485 124155
122265 139815
122368 123152
141664 153176
142310 168730
171856 176336
176272 180848
185368 203432
196724 202444
280540 365084
308620 389924
319550 430402
356408 399592
437456 455344
469028 486178
503056 514736
522405 525915
600392 669688
609928 686072
624184 691256
635624 712216
643336 652664
667964 783556
726104 796696
802725 863835
879712 901424
898216 980984
代码中的f就是7#的f
Amaris 发表于 2016-9-14 11:25
与 完整数 相关的是 亲和数。
用计算机找出10^6(即100万以内的亲和数对)，看看时间花多少？
耗时24小时哈哈哈费计算机啊 Everley 发表于 2016-9-14 11:28
纠正一下，取奇数或者偶数行没有办法解决问题2
所以代码改为
得
学习了好方法
页:
[1]