升序素数和降序素数

无心人

如果一个素数的十进制表示每个数位写成一个序列，其值是升序排列，则叫升序素数
否则叫降序素数
比如
  1123 1129就是升序素数
  54311就是降序素数

现在求100000000内的该类素数

medie2005

范围太小了吧,10^8内单调的数只有C(10+8-1,8)=C(17,9)=24310个,只要列举出这24310个数,再判素就ok了,应该可以在1秒内出解.
把范围改大一点才比较有难度,比如,改成10^25.

无心人

那考虑$10^32$以内的吧

medie2005

10^32内的单调数个数有C(32+10-1,32)=C(41,9)=350343565个.
我这边产生这350343565个数对应的组合大概要7秒,再由这些组合得到单调数.
再用Miller-Rabin测试判素,这些过程加起来,大概要半个钟头.

也可以做一些筛选,比如可以快速判断某个组合对应的单调数是否被小素数2,3,5,7,11整除,这样的筛选,大概可以将用Miller-Rabin测试判素的数减少5倍左右.

另外,可以取前N个素因子的乘积S,再对某个单调数a以gcd(S,a)来筛选,以减少Miller-Rabin判素算法的调用次数.

还有一种基于hash技术的做法,把32位分成两半,这样每一半只有C916+10-1,16)=C(25,9)=312455,这样每个单调数可表示为:xi*10^16+xj (xj的首数字大于xi的末数字),不过我没想到什么好办法来筛选,如果解决了这个问题,应该可以较大幅度地提高速度.

无心人

可以考虑优化了

medie2005

写了个小程序算了一下10^20内的降序素数,有488441个候选(只做了一次Miller-Rabin测试).输出文件有9M多,发不上来.
发一个10^15内的降序候选素数表.

无心人

你用PRIMO证明下
程序在我的那个素性和分解的帖子里

哈哈
估计要耗费100小时证明这么多素数

medie2005

发一下代码:

1. #include   <windows.h>
   
2. #include   <stdio.h>
   
3. #include   "gmp.h"
   
4. 
   
5. FILE *pf=fopen("data.txt","w");
   
6. 
   
7. const   int  L=8;
   
8. int     count=3;
   
9. 
   
10. unsigned prime[]={
    
11.         7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,
    
12.                 73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,
    
13.                 179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,
    
14.                 283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,
    
15.                 419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,
    
16.                 547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,
    
17.                 661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,
    
18.                 811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
    
19.                 947,953,967,971,977,983,991,997
    
20. };
    
21. 
    
22. unsigned p_prime[]={
    
23.         7,11,31,41,43,53,61,71,73,83,97,211,311,331,
    
24.                 421,431,433,443,521,541,557,631,641,643,653,
    
25.                 661,733,743,751,761,773,811,821,853,863,877,
    
26.                 881,883,887,911,941,953,971,977,983,991,997
    
27. };
    
28. 
    
29. mpz_t   Hash[10][L], Prime_Product, com_factor;
    
30. 
    
31. void    Init( ){
    
32.         int i, j;
    
33.         for( i=0; i<10; ++i ){
    
34.                 for( j=0; j<L; ++j )
    
35.                         mpz_init_set_ui( Hash[i][j], (i==0)? 0 : i );
    
36.         }
    
37.         for( i=1; i<10; ++i ){
    
38.                 for( j=1; j<L; ++j )
    
39.                         mpz_mul_ui( Hash[i][j], Hash[i][j-1], 10 );
    
40.         }       
    
41.         mpz_init_set_ui( Prime_Product, 1 );
    
42.         mpz_init( com_factor );
    
43.         for( i=0; i<sizeof(prime)/sizeof(prime[0]); ++i )
    
44.                 mpz_mul_ui( Prime_Product, Prime_Product, prime[i] );
    
45. }
    
46. 
    
47. void    Clear( ){
    
48.         mpz_clear( Prime_Product );
    
49.         mpz_clear( com_factor );
    
50.         for( int i=0; i<10; ++i ){
    
51.                 for( int j=0; j<L; ++j )
    
52.                         mpz_clear( Hash[i][j] );
    
53.         }
    
54.         fclose( pf );
    
55. }
    
56. 
    
57. void    Specail_Print(  ){
    
58.         fprintf(pf,"2\n3\n5\n");
    
59.         for( int k=0; k<sizeof(p_prime)/sizeof(p_prime[0]); ++k )
    
60.                 fprintf(pf,"%u\n",p_prime[k]), ++count;
    
61. }
    
62. 
    
63. void     combination(   int   n,   int   m   ){
    
64.         mpz_t num;
    
65.         mpz_init_set_ui( num, 0 );
    
66.         unsigned  *c=new unsigned[m+2], j, dig_sum=m;
    
67.         for( j=1; j<=m; ++j )  c[j]=j-1, mpz_add( num, num, Hash[1][j-1] );
    
68.         c[m+1]=n;       
    
69. R2:
    
70.         if( c[1]%2==0 && c[1]+1!=5 && dig_sum%3!=0 ){
    
71.                 mpz_gcd( com_factor, num, Prime_Product );               
    
72.                 if( mpz_cmp_ui( com_factor, 1 )==0 && mpz_probab_prime_p( num, 1 ) ){                       
    
73.                         mpz_out_str( pf, 10, num );
    
74.                         fprintf(pf,"\n");
    
75.                         ++count;
    
76.                 }
    
77.         }       
    
78.         if( m&1 )
    
79.                 if( c[1]+1<c[2] ){
    
80.                         ++c[1];
    
81.                         mpz_add_ui( num, num, 1 );
    
82.                         ++dig_sum;
    
83.                         goto R2;
    
84.                 }
    
85.                 else{
    
86.                         j=2;  goto R4;
    
87.                 }
    
88.                 else
    
89.                         if( c[1]>0 ){                       
    
90.                                 --c[1];
    
91.                                 mpz_sub_ui( num, num, 1 );
    
92.                                 --dig_sum;
    
93.                                 goto R2;
    
94.                         }
    
95.                         else{
    
96.                                 j=2;  goto R5;
    
97.                         };
    
98. R4:  if( c[j]>=j ){
    
99.                  if( j<=m )
    
100.                          mpz_sub( num, num, Hash[c[j]-j+2][j-1] ), dig_sum-=c[j]-j+2;                         
     
101.                  if( j-1<=m )
     
102.                          mpz_sub( num, num, Hash[c[j-1]-j+3][j-2] ), dig_sum-=c[j-1]-j+3;
     
103.          c[j]=c[j-1];
     
104.                  c[j-1]=j-2;
     
105.                  if( j<=m )
     
106.                          mpz_add( num, num, Hash[c[j]-j+2][j-1] ), dig_sum+=c[j]-j+2;
     
107.                  if( j-1<=m )
     
108.                          mpz_add( num, num, Hash[c[j-1]-j+3][j-2] ), dig_sum+=c[j-1]-j+3;
     
109.                  goto R2;
     
110.      }
     
111.      else
     
112.                  ++j;
     
113. R5:  if( c[j]+1<c[j+1] ){
     
114.                  if( j<=m )
     
115.                          mpz_sub( num, num, Hash[c[j]-j+2][j-1] ), dig_sum-=c[j]-j+2;
     
116.                  if( j-1<=m )
     
117.                          mpz_sub( num, num, Hash[c[j-1]-j+3][j-2] ), dig_sum-=c[j-1]-j+3;
     
118.          c[j-1]=c[j];
     
119.                  c[j]=c[j]+1;
     
120.                  if( j<=m )
     
121.                          mpz_add( num, num, Hash[c[j]-j+2][j-1] ), dig_sum+=c[j]-j+2;
     
122.                  if( j-1<=m )
     
123.                          mpz_add( num, num, Hash[c[j-1]-j+3][j-2] ), dig_sum+=c[j-1]-j+3;
     
124.                  goto R2;
     
125.      }
     
126.      else{
     
127.          ++j;
     
128.          if( j<=m )    goto R4;
     
129.      }
     
130.      delete []c;
     
131.          mpz_clear( num );
     
132. }
     
133. 
     
134. int   main(int   argc,   char   *argv[])
     
135. {       
     
136.         Init( );
     
137.         Specail_Print( );       
     
138.         int time=GetTickCount();
     
139.         for( int i=4; i<=L; ++i )
     
140.                 combination(   i+8,  i  );
     
141.         Clear( );
     
142.         printf("total : %d\n",count);
     
143.         printf("%d ms\n",GetTickCount()-time);
     
144.         return   0;
     
145. }

复制代码

medie2005

另外,如果要验证我在6#给的数据是否全是素数,可以采用特殊底的Miller-Rabin测试,底只需要取成2,3,7,61,24251这几个素数就可以了.
我已经用上面的方法对6#的数据做了全面的验算,表中的数据全部都是素数.

无心人

效果不错
程序啰嗦

哈