声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2314|回复: 11

[编程技巧] 关于完全数的matlab程序设计

[复制链接]
发表于 2016-9-14 11:20 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
来自:百度百科简介  完全数(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. a=1:10000;
  2. 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))))
复制代码


回复
分享到:

使用道具 举报

发表于 2016-9-14 11:21 | 显示全部楼层

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

    28

   496

        8128

Elapsed time is 7.991244 seconds.       我的电脑太慢了。。 .循环液尽量少用,费时!~
发表于 2016-9-14 11:22 | 显示全部楼层
解法二,用的是完整数的公式。
受他启发,直接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)
 楼主| 发表于 2016-9-14 11:23 | 显示全部楼层
上次在想如何去掉数组中的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))));
复制代码


发表于 2016-9-14 11:25 | 显示全部楼层
与 完整数 相关的是 亲和数。

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

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


点评

耗时24小时 哈哈哈 费计算机啊  详情 回复 发表于 2016-9-20 09:08
 楼主| 发表于 2016-9-14 11:25 | 显示全部楼层
1.请看我关于完全数的第三种做法,就是用向量操作计算约数和的另一种方法,就是用factor实现的,基于完全数的程序,加几行代码就可以实现亲和数
2.只需要取所得数的奇数行或者偶数行
因此我的代码如下

  1. clear;clc;close all
  2. clear
  3. a=1:10000;
  4. 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';
  5. e=ones(max(d),1);
  6. e(a)=d;e(1)=1;
  7. f=[e(e(e(a))==a'),e(e(e(e(a))==a'))];
  8. f(f(:,1)==f(:,2),:)=[];
  9. g=f(2:2:end,:)
复制代码

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

复制代码


 楼主| 发表于 2016-9-14 11:26 | 显示全部楼层
  1. g =

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

  42. Elapsed time is 1217.180373 seconds.
复制代码
这是在我的电脑上运行1000000内的结果
发表于 2016-9-14 11:27 | 显示全部楼层
Everley 发表于 2016-9-14 11:25
1.请看我关于完全数的第三种做法,就是用向量操作计算约数和的另一种方法,就是用factor实现的,基于完全数 ...

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


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

如7#的25~28行,等

87633, 69615
88730, 79750
100485, 124155
122265, 139815
 楼主| 发表于 2016-9-14 11:28 | 显示全部楼层
纠正一下,取奇数或者偶数行没有办法解决问题2
所以代码改为

  1. clear;clc;close all
  2. clear
  3. a=1:10000;
  4. 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';
  5. e=ones(max(d),1);
  6. e(a)=d;e(1)=1;
  7. f=[e(e(e(a))==a'),e(e(e(e(a))==a'))];
  8. f(f(:,1)==f(:,2),:)=[];
  9. h=sum(f,2);
  10. [n, m] = unique(h, 'first');
  11. p=f(sort(m),:)
复制代码

  1. p =

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


点评

学习了 好方法  详情 回复 发表于 2016-9-20 09:09
 楼主| 发表于 2016-9-14 11:29 | 显示全部楼层
一种更好的处理办法如下

  1. >> q=unique(sort(f,2),'rows')
  2. q =
  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.        63020       76084
  12.        66928       66992
  13.        67095       71145
  14.        69615       87633
  15.        79750       88730
  16.       100485      124155
  17.       122265      139815
  18.       122368      123152
  19.       141664      153176
  20.       142310      168730
  21.       171856      176336
  22.       176272      180848
  23.       185368      203432
  24.       196724      202444
  25.       280540      365084
  26.       308620      389924
  27.       319550      430402
  28.       356408      399592
  29.       437456      455344
  30.       469028      486178
  31.       503056      514736
  32.       522405      525915
  33.       600392      669688
  34.       609928      686072
  35.       624184      691256
  36.       635624      712216
  37.       643336      652664
  38.       667964      783556
  39.       726104      796696
  40.       802725      863835
  41.       879712      901424
  42.       898216      980984
复制代码
代码中的f就是7#的f
发表于 2016-9-20 09:08 | 显示全部楼层
Amaris 发表于 2016-9-14 11:25
与 完整数 相关的是 亲和数。

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

耗时24小时  哈哈哈  费计算机啊
发表于 2016-9-20 09:09 | 显示全部楼层
Everley 发表于 2016-9-14 11:28
纠正一下,取奇数或者偶数行没有办法解决问题2
所以代码改为

学习了  好方法
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-23 06:19 , Processed in 0.065555 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表