[讨论]:查找所有“循环质数”

bainhome

cody原帖请见   Problem 1095. Circular Primes 。
实际这仍然是Project Euler中的问题(Pro.35:Circular Prime)，移植到MATLAB来而已。
原题内容如下：

1. The number, 197, is called a circular prime because all rotations of the digits: 197, 971, and 719, are themselves prime.
   
2. 
   
3. There are thirteen such primes below 100: 2, 3, 5, 7, 11, 13, 17, 31, 37, 71, 73, 79, and 97.
   
4. 
   
5. Given a number x, write a MATLAB script that will tell you the number of circular primes less than or equal to x as well as a sorted list of what the circular prime numbers are.
   

复制代码

意思给定一个整数N，找出所有小于等于该整数的循环质数序列，所谓循环质数是指这个质数所有位数上的数字循环旋转(注意：不是所有位数上数字的排列或者组合！)都是质数，例如：

1. x=[197,719,971]

复制代码

上述就是满足条件的三个循环质数，但如：

1. x_i=[791,179]

复制代码

或会因不符合循环旋转条件、或会因不是质数而不在序列内。
这个问题比较有意思的地方我个人感觉有两个：
1.如何做数字换位，大家看看有多少不同思路？
2.如何在保证效率的条件下减小size
这次诸位把自己的代码贴出来以我的电脑为准测试运行时间，祁工算一下size(暂时不要公布这个size的计算办法，虽然很多人已经知道，不过还是越少越好，这真心不是个好东西)，哥儿几个休息好了没有？开工了...
最终讨论个size和运行时间的权重计算办法，按该办法得到的最优代码额外给个技术分。大家看看是否可行？

liuyalong008

我还是老老实实贴个长的吧

1. function [how_many what_numbers]=circular_prime(x)
   
2. % x=1000
   
3. what_numbers=[2;3;5];
   
4. m=reshape(sprintf('%10d',(7:x)'),10,[])';
   
5. r=str2num(m((sum(m=='1'|m=='3'|m=='7'|m=='9'|m==' ',2)==10),:));
   
6. f=r(isprime(r));
   
7. for k=1:length(f)
   
8.     s=num2str(f(k))-'0';
   
9.     len=numel(s);
   
10.     t=repmat(1:len,len,1);
    
11.     t1=t+rot90(t)-1;
    
12.     t1(t1>len)=t1(t1>len)-len;
    
13.     a = str2num(char(s(t1)+'0'));
    
14.     what_numbers=[what_numbers;a(isprime(a)&all(isprime(a)))];
    
15. end
    
16. what_numbers=unique(what_numbers);
    
17. what_numbers(what_numbers>x)=[];
    
18. how_many=numel(what_numbers);

复制代码

期待看到精彩的

liuyalong008

再贴一个去掉循环后的结果：
此次速度明显降低，估计size和上个相当

1. function[how_many, what_numbers]=circular_prime_1(x)
   
2. % x=54321
   
3. num=(7:x)';
   
4. m=reshape(sprintf('%5d',num),5,[])';
   
5. index=cell2mat(arrayfun(@(i)m(:,i)==49|m(:,i)==51|m(:,i)==55|m(:,i)==57|m(:,i)==32,1:5,'Uni',0));
   
6. m=m(isprime(num)&sum(index,2)==5,:);
   
7. t1=cell2mat(arrayfun(@(i)circshift(1:5,[0 ,i]),(1:5)','Uni',0));
   
8. s=char(arrayfun(@(i)m(:,t1(i,:)),1:5,'Uni',0));
   
9. data=str2num(char(regexprep(num2cell(s,2),' ','')));
   
10. r=isprime(data);
    
11. d=reshape(data,[],5);
    
12. what_numbers=[2; 3; 5; d(sum(reshape(r,[],5),2)==5,5)];
    
13. how_many=numel(what_numbers);
    

复制代码

nwcwww

1. function [how_many ans]=circular_prime(x)
   
2. str2double(regexp(num2str(1:2:x),'\<[1379]+\>','match'));
   
3. ans(isprime(ans));
   
4. strtrim(cellstr(num2str(ans')));
   
5. anon = @(y) transpose(str2num(y(mod(bsxfun(@plus,(0:length(y)-1)',0:length(y)-1),length(y))+1)));
   
6. FullList = cellfun(anon, ans, 'uni', 0);
   
7. IfCirP = cellfun(@(a) all(isprime(a)), FullList);
   
8. unique([2 5 FullList{IfCirP}]);
   
9. ans(ans<=x);
   
10. how_many = length(ans);
    
11. end

复制代码

size = 118

如果考虑size的问题的话，大概能减到85左右吧。再删减可能得换思路了。

qibbxxt

nwcwww 发表于 2013-6-8 08:04
size = 118

如果考虑size的问题的话，大概能减到85左右吧。再删减可能得换思路了。 ...

ncw 老兄的思维很活跃啊
在此基础上改一个

1. function [i,ans]=circular_prime(x)
   
2. a = regexp(num2str(primes(x)),'\<[1379]+\>','match');
   
3. cellfun(@(x)all(isprime(gallery('circul',x-'0')*10.^(length(x)-1:-1:0)')),a(3:end));
   
4. unique([primes(10),str2num(char(a([true true ans])))']);
   
5. i = numel(ans);
   
6. end

复制代码

rocwoods

我贴个size198的，比较适合较大规模x的代码(x = 9e6)，避免字符串转来转去，过滤掉部分容易判断的非素数，size优化可以在这个基础上进行：

1. function [how_many,num] = myCircuPrime(x)
   
2. t = [1 3 7 9];
   
3. num = [2 3 5 7];
   
4. for k = 2:floor(log10(x))+1
   
5.     allnum = t(fullfact(4*ones(1,k)));
   
6.     allnumD = sum(bsxfun(@times,allnum,10.^(k-1:-1:0)),2);
   
7.     indp1 = find(all(bsxfun(@mod,allnumD,[3 7 11 13 17]),2));
   
8.     indp2 = find(isprime(allnumD(indp1)));
   
9.     allnum = allnum(indp1(indp2),:);
   
10.     ind = gallery('circul',1:k);
    
11.     for kk = 1:length(allnum)
    
12.         tmp = allnum(kk,:);
    
13.         if all(isprime(sum(bsxfun(@times,tmp(ind),10.^(k-1:-1:0)),2)))
    
14.             num(end+1) = allnumD(indp1(indp2(kk)));
    
15.         end
    
16.     end
    
17. end
    
18. num = sort([11 13 17,num(num<=x)]);
    
19. how_many = length(num);

复制代码

bainhome

大家发现没有:吴某人代码非常喜欢用最低层的基本命令，这使得代码执行效率、我是说执行效率，非常高——我现在最想实现的就是这种意图。

lin2009

各位八仙过海各显神通，但是没有注释的代码看起来很费劲，特别是有的函数不是Matlab的基础函数（如fullfact）
大家都忽略了这是个输出个数有限的问题。n位数的循环质数如下：
bit1 = 2,3,5,7,
bit2 = 11,13,17,31,37,71,73,79,97,
bit3 = 113,131,197,199,311,337,373,719,733,919,971,991,
bit4 = 1193,1931,3119,3779,7793,7937,9311,9377,
bit5 = 11939,19391,19937,37199,39119,71993,91193,93719,93911,99371,
bit6 = 193939,199933,319993,331999,391939,393919,919393,933199,939193,939391,993319,999331
bit7 = 7位数：无
bit8 = 8位数：无
bit9 = 9位及更高位数的没有涉及到，猜测循环质数应该不存在。
因为个数有限（至少在9位数以下），因此从解决问题的实用角度可以先找出某一范围内所有的循环质数再查表：
[code]
function [nums,output] = myprimes(x)
Circle_primes = [...
    2,3,5,7,...
    11,13,17,31,37,71,73,79,97,...
    113,131,197,199,311,337,373,719,733,919,971,991,...
    1193,1931,3119,3779,7793,7937,9311,9377,...
    11939,19391,19937,37199,39119,71993,91193,93719,93911,99371,...
    193939,199933,319993,331999,391939,393919,919393,933199,939193,939391,993319,999331];
output = Circle_primes( x >= Circle_primes ); % 输出比x小的循环质数。
nums = length(output);
[\code]
从决解问题的角度看，这个“算法”最快，Just for fun。
(Given a number x, write a MATLAB script that will tell you the number of circular primes less than or equal to x
as well as a sorted list of what the circular prime numbers are.)

这个问题（如上）应该描述为在某个范围内寻找这些循环质数的的最快算法。（其实大家都是这样做的）
很显然要判断质数，程序绕不开isprime这个耗时的函数(不大可能自己搞个什么算法)，因此尽量减少需要判断质数的数量。
最快的是用[1,3,7,9]这几个数字（可以直接是double等数值类型，不必通过char转换）进行排列组合。（显然其它数字不可能构成循环质数）
再在这基础上删除能被3整除的数（或者再根据什么规律删除其他非质数，尽量不用isprime来判断），
减少需要调用isprime的个数。然后再进行判断。

另外，程序运行时间除了与机器配置等有关外，还与程序本身的具体参数有关，及myprimes(x)中参数x的值。
显然x = 100时， x = 1e6时 和 x = 1e8 时的运行时间是有很大的差别的。
没有标明这些，单独列出的运行时间是没有太大的参考价值的，即使是同一台机子。

bainhome

如下点评和总结全部来自祁工，这货这两天比较忙，我帮他发掉吧。(所以别给我加分，我在上面帖子里给他加一下即可)。
他一贯的风格：一句废话没有，朴实、心无旁骛直冲主题。
一、循环质数这道题目的基本思路：

思路1

1. 找出x之前的所有质数

2. 逐个判断质数循环是否满足条件

3. 去掉不满足条件的或者加入满足条件者

思路2

1. 观察规律，最后的答案出现的数字都是1,3,7,9；因此先生成1,3,7,9的各种组合

2. 逐个判断是否满足条件

3. 去掉不满足条件的或者加入满足条件者

对比：前者用primes函数节省了人工观察规律的部分时间，将任务交给计算机；后者是从算法上改进，后者的思路还是很可取的，有一个规律，想出第二种思路的都是中国人。

二、实现细节

1. 找出x之前的所有质数,用primes函数来完成。

2. 提取数字，常见的有两种思路

         2.1将数字转化为字符串，利用字符串的ACSII码来转化，这应该是属于一个技巧

1. int2str(x)
   
2. - '0'

复制代码

         2.2利用公式来求

1. mod(floor(x./10.^(n-1:-1:0)),10)

复制代码

3. 数字循环

         3.1利用循环

1. x= [x(2:end),x(1)];

复制代码

         3.2利用circshift

1. circshift(x,k)

复制代码

                   注意：x和k的关系

         3.3用测试矩阵

1. d=gallery('circul',x);

复制代码

         3.4其他函数，比如hankel等

         3.5构造，比如liu的rot90

4. 有向量转化为数字

         4.1常见

1. str2num

复制代码

         4.2用

1. polyval(x,10);

复制代码

         4.3用矩阵乘法，

1. x*10^(0:n-1)

复制代码

5. 判断是否为质数

         5.1isprime

         5.2向量判断all,any的使用

6. 删除重复元素函数unique

三、部分数组组成的数字的计算方法

1.使用排列组合的方法

         1.1使用fullfact函数

         如生成[1 3 7 9]组成的五位数

1. a = [1 3 7 9];
   
2. b = 10.^(0:4)*a(fullfact(numel(a)+zeros(1,5)))'

复制代码

         1.2使用combvec函数

1. b= 10.^(0:4)*combvec(a,a,a,a,a);

复制代码

2.使用字符串操作

         2.1用正则表达式代替

1. str2double(regexp(num2str(1:2:10^ans),'\<[1379]+\>','match'));%（nwc）
   
2. unique(str2num(regexprep(num2str(1:100),'[024568]','')))

复制代码

         2.2查找

1. r=str2num(m((sum(m=='1'|m=='3'|m=='7'|m=='9'|m==' ',2)==10),:));%(Liu)

复制代码

         这个‘’的技巧很巧妙

四 arrayfun,cellfun的使用

         这写函数的使用需要对矩阵的形式构造，尤其是这个问题，由于数字的位数不同，导致数组的统一操作有点麻烦，所以构造arrayfun，去掉for循环，就需要一些技巧

五leading

leadingsolution 来自cody中大家都崇拜的大神Alfonso Nieto-Castanon

1. 
   
2. function [ans,i]=circular_prime(x)
   
3. primes(x);
   
4. i=ans(arrayfun(@(n)all(arrayfun(@(k)isprime(str2num(circshift(num2str(n)',k)')),1:4)),ans));
   
5. numel(i);
   
6. end
   

复制代码

对了，大家讨论一下最优思路，给个技术分，没那么严格，本来就是玩儿，祁工qq上说权重方式不一定好，可以用帕累托集合度量，我是不会，交给你度量吧，哈哈。

bainhome

我只想说leading solution这种以size为评价标准的做法真心非常不好！这个帖子就是最佳例子。如下是leading solution的运行时间：

1. tic;[ans,i]=circular_prime_alfonso(54321),toc;
   
2. ans =
   
3.     38
   
4. i =
   
5.   Columns 1 through 12
   
6.            2           3           5           7          11          13          17          31          37          71          73          79
   
7.   Columns 13 through 24
   
8.           97         113         131         197         199         311         337         373         719         733         919         971
   
9.   Columns 25 through 36
   
10.          991        1193        1931        3119        3779        7793        7937        9311        9377       11939       19391       19937
    
11.   Columns 37 through 38
    
12.        37199       39119
    
13. Elapsed time is 4.599584 seconds.

复制代码

毋庸置疑，Alfonso Nieto-Castanon函数操控水准绝对属顶级大师，可在size至上的规则下，他都不能按问题驱动的思路构造真正“最有效率”的算法流程，遑论我们。
而这一点，无疑是cody中的一个无法弥补的遗憾。

bainhome

另外，刚才和祁工一起检讨了一下，的确如lin2009所说，我们大家代码不加注释并不是一种很好的习惯。让其他人很难阅读，今后改正。