数学研发网设为首页收藏本站
切换到窄版

 找回密码
 欢迎注册
数学研发论坛»论坛 › 【算法研发】 › 编程擂台 › 估计下找到10^100以上的一个四生素数的工作量 ...
1234567下一页
返回列表 发新帖
楼主: 无心人

[讨论] 估计下找到10^100以上的一个四生素数的工作量

[复制链接]
mathe
发表于 2008-2-14 16:57:50 | 显示全部楼层
奇怪，实际计算表明用TestPrimality(5)去预选，又慢了很多，郭老大能猜测到底是什么原因吗？感觉挺奇怪的
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

mathe
发表于 2008-2-14 17:29:15 | 显示全部楼层
还有我上面测试10次，平均每次得到4.4组结果。 根据前面楼主的公式，四生素数分布密度大概为$\frac{4.1}{log^4(x)}$ 由于程序里面使用一个$10^7$大小的数组来筛选数据，数组每项代表32个位置，每个位置包含30个数（之看模30余15的情况） 所以总共可以覆盖$10^7 xx 32 xx 30$个数据，所以平均应该有$10^7 xx 32 xx30xx\frac{4.1}{log^4(10^100)} = 14.0$ 这个偏差挺大的，也许楼主给出公式里面每组数看成4个，这样平均应该3.5个数，结果就比较相近了。
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

mathe
发表于 2008-2-14 18:23:26 | 显示全部楼层
发现我的代码的一处BUG. 下面代码中int s;应该改成long long s; 要不然在SQRTSMALLPRIMES>215时计算就有可能发生溢出，从而出现错误的计算，导致结果偏少。 修改以后，明显找到的结果要增加很多。 还有一点发现的结果是SQRTSMALLPRIMES的增加对速度影响不是很大，所以可以适当提高，我的机器上增加到4000后计算速度 可以有显著提高，总运算时间基本都在2分钟之内（筛选1分钟，后面测试1分钟）
for(k=0;k
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

mathe
发表于 2008-2-14 18:23:58 | 显示全部楼层
新的代码:
  1. // HugeCalcDemo.cpp : Defines the entry point for the console application.
  2. //
  4. // Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
  5. // Win32 Debug: Debug Multithreaded DLL
  6. // Win32 Release: Multithreaded DLL
  8. #include <time.h>
  9. #include <iostream>
  10. using namespace std;
  11. #include "../../../../HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
  12. #include "../../../../HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
  13. #include "../../../../HugeCalc_API/CppAPI/Include/HugeIntX.h" // 16进制系统
  15. #pragma message( "automatic link to ../../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
  16. #pragma comment( lib, "../../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
  18. // http://www.swox.com/gmp/#TRY
  19. // 1 + gcd(87324,78263148,7896) * (10^1989879887 mod 471!)
  21. #if 1
  22. #define _Test_HI
  23. #define integer CHugeInt
  24. #else
  25. #define _Test_HX
  26. #define integer CHugeIntX
  27. #endif
  29. #define SQRTSMALLPRIMES 4000
  30. #define SMALLPRIMES (SQRTSMALLPRIMES*SQRTSMALLPRIMES)
  31. int spcount;
  32. int sptable[SMALLPRIMES];
  34. #define LIMIT_SIZE 10000000
  35. #define RANGE (LIMIT_SIZE<<5)
  36. int mask[LIMIT_SIZE];
  37. #define SET_MASK(x) (mask[(x)>>5]|=1<<((x)&31))
  38. #define IS_SET(x) (mask[(x)>>5]&(1<<((x)&31)))
  40. void initsp()
  41. {
  42. int i,j;
  43. memset(sptable,-1,sizeof(sptable));
  44. sptable[0]=sptable[1]=0;
  45. for(i=2;i<SQRTSMALLPRIMES;i++){
  46. if(sptable[i]){
  47. for(j=i*i;j<SMALLPRIMES;j+=i){
  48. sptable[j]=0;
  49. }
  50. }
  51. }
  52. for(i=0,spcount=0;i<SMALLPRIMES;i++){
  53. if(sptable[i]){
  54. sptable[spcount++]=i;
  55. }
  56. }
  57. }
  60. int main(int argc, char* argv[])
  61. {
  62. // printf("Hello World!\n");
  64. cout << "Call " << HugeCalc::GetVer() << endl << endl;
  66. if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel())
  67. {
  68. cout << endl << "警告：您未通过 HugeCalc.dll 的许可认证！" \
  69. << endl << endl << "解决方案可选下列方案之一：" \
  70. << endl << " 一、请将本程序移动到“/CopyrightByGuoXianqiang/[../]”目录下运行；" \
  71. << endl << " 二、或请在 HugeCalc.chm 中进行注册（一劳永逸）。" \
  72. << endl << endl;
  73. system( "pause" );
  74. return (-1);
  75. }
  76. // 初始化
  77. HugeCalc::EnableTimer( TRUE );
  78. HugeCalc::ResetTimer();
  79. HugeCalc::SeedRandom(time(NULL));
  81. initsp();
  83. integer hugeStart, hugeTemp;
  84. #ifdef _Test_HI
  85. hugeStart.Random(100);
  86. #else
  87. hugeStart.Random(335);
  88. #endif
  89. int m=hugeStart%30;
  90. m=15-m;
  91. if(m<0)m+=15;
  92. hugeStart+=m;///After the normalization, hugeStart%30==15
  94. int filter_list[]={-4,-2,2,4};
  95. int i;
  96. for(i=3;i<spcount;i++){//skip prime 2,3,5
  97. int p=sptable[i];
  98. m=hugeStart%p;
  99. int v=HugeCalc::InvertMod(30,p);
  100. int k;
  101. for(k=0;k<sizeof(filter_list)/sizeof(filter_list[0]);k++){
  102. long long s=-m+filter_list[k];
  103. while(s<0)s+=p;
  104. s*=v;
  105. s%=p;//hugeStart+s*30==filter_list[k] (mod p)
  106. int h;
  107. for(h=(int)s;h<RANGE;h+=p){
  108. SET_MASK(h);///Filter hugeStart+h*30
  109. }
  110. }
  112. }
  114. CONST UINT32 u32Timer_Filter = HugeCalc::GetTimer( TIMER_UNIT_ms );
  115. cout << "Filter cost "<<u32Timer_Filter<<"ms"<<endl;
  116. HugeCalc::ResetTimer();
  117. hugeStart-=4;
  119. int tcount=0,gcount=0;
  121. for(i=0;i<RANGE;i++){
  122. if(!IS_SET(i)){///hugeStart+i*30 may be a candidate
  123. tcount++;
  124. #if 0
  125. if(((hugeTemp=hugeStart)).TestPrimality(5)!=COMPOSITE&&
  126. (hugeTemp+=2).TestPrimality(5)!=COMPOSITE&&
  127. (hugeTemp+=4).TestPrimality(5)!=COMPOSITE&&
  128. (hugeTemp+=2).TestPrimality(5)!=COMPOSITE){
  129. #endif
  130. if(((hugeTemp=hugeStart)).IsPrime()&&
  131. (hugeTemp+=2).IsPrime()&&
  132. (hugeTemp+=4).IsPrime()&&
  133. (hugeTemp+=2).IsPrime()){
  134. ///Output the candidate;
  135. CONST LPCTSTR lpShow = hugeStart.GetStr( FS_DEFAULT );
  136. gcount++;
  137. cout<<"Found n,n+2,n+6,n+8 where n is "<<lpShow <<endl;
  138. }
  139. #if 0
  140. }
  141. #endif
  143. }
  144. hugeStart+=30;
  145. }
  148. // 记录计算耗时
  149. CONST UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_ms );
  152. // 输出
  153. cout << "Total took " << u32Timer_Calc << " ms" << endl;
  154. cout << "Total test "<<tcount<<" Numbers, get "<<gcount<<" numbers"<<endl;
  156. system( "pause" );
  158. return 0;
  159. }
复制代码
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

gxqcn
发表于 2008-2-14 18:30:43 | 显示全部楼层
原帖由 mathe 于 2008-2-14 16:57 发表 奇怪，实际计算表明用TestPrimality(5)去预选，又慢了很多，郭老大能猜测到底是什么原因吗？感觉挺奇怪的
Rabin-Miller Strong Pseudoprime Test 误判概率约为$1/4^N$或更低，其中N为随机测试次数。 TestPrimality(1) 误判概率约为$1/4^1 = 25%$，显得太弱了； TestPrimality(5) 误判概率约为$1/4^5 = 0.0976%$，已经很强了，许多软件基本“认为”其是个可信的“素数”了。 所以，要适当掌握度，也许 TestPrimality(3) 是个不错的选择？
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

mathe
发表于 2008-2-14 18:35:52 | 显示全部楼层
修改BUG后的程序运行数次分别得到数目为 15,0,16,18,16 所以数目同估计公式还是基本匹配的
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

mathe
发表于 2008-2-15 10:17:03 | 显示全部楼层
Faint:
  1. if(m<0)m+=15;
复制代码
shoulbe be changed to=>
  1. if(m<0)m+=30;
复制代码
That's why there're so much zeros.
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

gxqcn
发表于 2008-2-15 10:37:45 | 显示全部楼层
我将19#重新编辑了下。 将 mathe 代码中的 initsp() 替换成 "HugeCalc::GetPrimeList( sptable, spcount, 0, SMALLPRIMES );" 并将所有 int 型变量替换成无符号型；以及其它一些小改动。。。
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

mathe
发表于 2008-2-15 10:49:04 | 显示全部楼层
Another small optimization opportunity: change loop
  1. for(i=0;i<RANGE;i++){
  2. if(!IS_SET(i)){//hugeStart+i*30 may be a candidate
  3. ...
  4. }
  5. }
复制代码
to
  1. for(i=0;i<LIMIT_SIZE;i++){
  2. if(!mask[i]){
  3. hugeStart+=30*32;
  4. }else{
  5. for(t=0;t<32;t++){
  6. if(mask[i]&(1<<t)){
  7. .....
  8. }
  9. hugeStart+=30;
  10. }
  11. }
  12. }
复制代码
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

gxqcn
发表于 2008-2-15 16:34:34 | 显示全部楼层
因为我们是从 p --> p+2 --> p+6 --> p+8 这个顺序逐步测试素性的，应从“粗放”到“精细”过渡， 所以排在前面的数应以淘汰效率为主；排在后面的数应以准确度为主，比方说最后一个数就应直接调用一次 IsPrime() 即可。 我将原来的代码：
  1. if( (hugeTest=hugeStart).IsPrime()&&
  2. (hugeTest+=2).IsPrime()&&
  3. (hugeTest+=4).IsPrime()&&
  4. (hugeTest+=2).IsPrime()){
  5. //Output the candidate;
复制代码
替换成新的（已在19#楼修正）：
  1. if( (hugeTest=hugeStart).TestPrimality(1) != COMPOSITE &&
  2. (hugeTest+=2).TestPrimality(2) != COMPOSITE &&
  3. (hugeTest+=4).TestPrimality(3) != COMPOSITE &&
  4. (hugeTest+=2).IsPrime() &&
  5. (hugeTest-=2).IsPrime() &&
  6. (hugeTest-=4).IsPrime() &&
  7. (hugeTest-=2).IsPrime()){
  8. //Output the candidate;
复制代码
搜索$10^100$后的连续13组“四生素数”，时间由先前的1'22.174''降为1’13.607''，提速约11.64% 补记：19:00 左右，再次对代码优化，去掉了临时变量“hugeTest”，节省了一次大数拷贝动作。但却用了一个C/C++少有的语法——在if逻辑语句中大量使用逗号“,”运算操作符。经测试，可在新标准下再提速18.31%
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

1234567下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-11-22 00:16 , Processed in 0.030960 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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