找回密码
 欢迎注册
查看: 23262|回复: 10

[擂台] 升序素数和降序素数

[复制链接]
发表于 2008-4-23 18:43:46 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有账号?欢迎注册

×
如果一个素数的十进制表示每个数位写成一个序列,其值是升序排列,则叫升序素数 否则叫降序素数 比如 1123 1129就是升序素数 54311就是降序素数 现在求100000000内的该类素数
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-4-23 21:25:28 | 显示全部楼层
范围太小了吧,10^8内单调的数只有C(10+8-1,8)=C(17,9)=24310个,只要列举出这24310个数,再判素就ok了,应该可以在1秒内出解. 把范围改大一点才比较有难度,比如,改成10^25.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-4-23 21:31:15 | 显示全部楼层
那考虑$10^32$以内的吧
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-4-24 10:38:16 | 显示全部楼层
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的末数字),不过我没想到什么好办法来筛选,如果解决了这个问题,应该可以较大幅度地提高速度.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-4-24 10:48:23 | 显示全部楼层
可以考虑优化了
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-4-24 21:39:26 | 显示全部楼层
写了个小程序算了一下10^20内的降序素数,有488441个候选(只做了一次Miller-Rabin测试).输出文件有9M多,发不上来. 发一个10^15内的降序候选素数表.

10^15内降序候选素数表.rar

183.89 KB, 下载次数: 8, 下载积分: 金币 -1 枚, 经验 1 点, 下载 1 次

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-4-24 22:11:06 | 显示全部楼层
你用PRIMO证明下 程序在我的那个素性和分解的帖子里 哈哈 估计要耗费100小时证明这么多素数
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-4-26 16:28:13 | 显示全部楼层
发一下代码:
  1. #include <windows.h>
  2. #include <stdio.h>
  3. #include "gmp.h"
  4. FILE *pf=fopen("data.txt","w");
  5. const int L=8;
  6. int count=3;
  7. unsigned prime[]={
  8. 7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,
  9. 73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,
  10. 179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,
  11. 283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,
  12. 419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,
  13. 547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,
  14. 661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,
  15. 811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
  16. 947,953,967,971,977,983,991,997
  17. };
  18. unsigned p_prime[]={
  19. 7,11,31,41,43,53,61,71,73,83,97,211,311,331,
  20. 421,431,433,443,521,541,557,631,641,643,653,
  21. 661,733,743,751,761,773,811,821,853,863,877,
  22. 881,883,887,911,941,953,971,977,983,991,997
  23. };
  24. mpz_t Hash[10][L], Prime_Product, com_factor;
  25. void Init( ){
  26. int i, j;
  27. for( i=0; i<10; ++i ){
  28. for( j=0; j<L; ++j )
  29. mpz_init_set_ui( Hash[i][j], (i==0)? 0 : i );
  30. }
  31. for( i=1; i<10; ++i ){
  32. for( j=1; j<L; ++j )
  33. mpz_mul_ui( Hash[i][j], Hash[i][j-1], 10 );
  34. }
  35. mpz_init_set_ui( Prime_Product, 1 );
  36. mpz_init( com_factor );
  37. for( i=0; i<sizeof(prime)/sizeof(prime[0]); ++i )
  38. mpz_mul_ui( Prime_Product, Prime_Product, prime[i] );
  39. }
  40. void Clear( ){
  41. mpz_clear( Prime_Product );
  42. mpz_clear( com_factor );
  43. for( int i=0; i<10; ++i ){
  44. for( int j=0; j<L; ++j )
  45. mpz_clear( Hash[i][j] );
  46. }
  47. fclose( pf );
  48. }
  49. void Specail_Print( ){
  50. fprintf(pf,"2\n3\n5\n");
  51. for( int k=0; k<sizeof(p_prime)/sizeof(p_prime[0]); ++k )
  52. fprintf(pf,"%u\n",p_prime[k]), ++count;
  53. }
  54. void combination( int n, int m ){
  55. mpz_t num;
  56. mpz_init_set_ui( num, 0 );
  57. unsigned *c=new unsigned[m+2], j, dig_sum=m;
  58. for( j=1; j<=m; ++j ) c[j]=j-1, mpz_add( num, num, Hash[1][j-1] );
  59. c[m+1]=n;
  60. R2:
  61. if( c[1]%2==0 && c[1]+1!=5 && dig_sum%3!=0 ){
  62. mpz_gcd( com_factor, num, Prime_Product );
  63. if( mpz_cmp_ui( com_factor, 1 )==0 && mpz_probab_prime_p( num, 1 ) ){
  64. mpz_out_str( pf, 10, num );
  65. fprintf(pf,"\n");
  66. ++count;
  67. }
  68. }
  69. if( m&1 )
  70. if( c[1]+1<c[2] ){
  71. ++c[1];
  72. mpz_add_ui( num, num, 1 );
  73. ++dig_sum;
  74. goto R2;
  75. }
  76. else{
  77. j=2; goto R4;
  78. }
  79. else
  80. if( c[1]>0 ){
  81. --c[1];
  82. mpz_sub_ui( num, num, 1 );
  83. --dig_sum;
  84. goto R2;
  85. }
  86. else{
  87. j=2; goto R5;
  88. };
  89. R4: if( c[j]>=j ){
  90. if( j<=m )
  91. mpz_sub( num, num, Hash[c[j]-j+2][j-1] ), dig_sum-=c[j]-j+2;
  92. if( j-1<=m )
  93. mpz_sub( num, num, Hash[c[j-1]-j+3][j-2] ), dig_sum-=c[j-1]-j+3;
  94. c[j]=c[j-1];
  95. c[j-1]=j-2;
  96. if( j<=m )
  97. mpz_add( num, num, Hash[c[j]-j+2][j-1] ), dig_sum+=c[j]-j+2;
  98. if( j-1<=m )
  99. mpz_add( num, num, Hash[c[j-1]-j+3][j-2] ), dig_sum+=c[j-1]-j+3;
  100. goto R2;
  101. }
  102. else
  103. ++j;
  104. R5: if( c[j]+1<c[j+1] ){
  105. if( j<=m )
  106. mpz_sub( num, num, Hash[c[j]-j+2][j-1] ), dig_sum-=c[j]-j+2;
  107. if( j-1<=m )
  108. mpz_sub( num, num, Hash[c[j-1]-j+3][j-2] ), dig_sum-=c[j-1]-j+3;
  109. c[j-1]=c[j];
  110. c[j]=c[j]+1;
  111. if( j<=m )
  112. mpz_add( num, num, Hash[c[j]-j+2][j-1] ), dig_sum+=c[j]-j+2;
  113. if( j-1<=m )
  114. mpz_add( num, num, Hash[c[j-1]-j+3][j-2] ), dig_sum+=c[j-1]-j+3;
  115. goto R2;
  116. }
  117. else{
  118. ++j;
  119. if( j<=m ) goto R4;
  120. }
  121. delete []c;
  122. mpz_clear( num );
  123. }
  124. int main(int argc, char *argv[])
  125. {
  126. Init( );
  127. Specail_Print( );
  128. int time=GetTickCount();
  129. for( int i=4; i<=L; ++i )
  130. combination( i+8, i );
  131. Clear( );
  132. printf("total : %d\n",count);
  133. printf("%d ms\n",GetTickCount()-time);
  134. return 0;
  135. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-4-26 16:50:18 | 显示全部楼层
另外,如果要验证我在6#给的数据是否全是素数,可以采用特殊底的Miller-Rabin测试,底只需要取成2,3,7,61,24251这几个素数就可以了. 我已经用上面的方法对6#的数据做了全面的验算,表中的数据全部都是素数.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-4-26 17:13:09 | 显示全部楼层
效果不错 程序啰嗦 哈
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-11-23 17:37 , Processed in 0.026896 second(s), 19 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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