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

mathe

奇怪，实际计算表明用TestPrimality(5)去预选，又慢了很多，郭老大能猜测到底是什么原因吗？感觉挺奇怪的

mathe

还有我上面测试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

发现我的代码的一处BUG.
下面代码中int s;应该改成long long s;
要不然在SQRTSMALLPRIMES>215时计算就有可能发生溢出，从而出现错误的计算，导致结果偏少。
修改以后，明显找到的结果要增加很多。
还有一点发现的结果是SQRTSMALLPRIMES的增加对速度影响不是很大，所以可以适当提高，我的机器上增加到4000后计算速度
可以有显著提高，总运算时间基本都在2分钟之内（筛选1分钟，后面测试1分钟）

        for(k=0;k<sizeof(filter_list)/sizeof(filter_list[0]);k++){
            int s=-m+filter_list[k];
            while(s<0)s+=p;
            s*=v;
            s%=p;//hugeStart+s*30==filter_list[k] (mod p)
            int h;
            for(h=s;h<RANGE;h+=p){
                SET_MASK(h);///Filter hugeStart+h*30
            }
        }

mathe

新的代码:

1. // HugeCalcDemo.cpp : Defines the entry point for the console application.
   
2. //
   
3. 
   
4. //        Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
   
5. //                Win32 Debug:        Debug Multithreaded DLL
   
6. //                Win32 Release:        Multithreaded DLL
   
7. 
   
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进制系统
    
14. 
    
15. #pragma message( "automatic link to ../../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
    
16. #pragma comment( lib, "../../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
    
17. 
    
18. //        http://www.swox.com/gmp/#TRY
    
19. //        1 + gcd(87324,78263148,7896) * (10^1989879887 mod 471!)
    
20. 
    
21. #if 1
    
22.         #define _Test_HI
    
23.         #define integer CHugeInt
    
24. #else
    
25.         #define _Test_HX
    
26.         #define integer CHugeIntX
    
27. #endif
    
28. 
    
29. #define SQRTSMALLPRIMES 4000
    
30. #define SMALLPRIMES (SQRTSMALLPRIMES*SQRTSMALLPRIMES)
    
31. int spcount;
    
32. int sptable[SMALLPRIMES];
    
33. 
    
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)))
    
39. 
    
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. }
    
58. 
    
59. 
    
60. int main(int argc, char* argv[])
    
61. {
    
62. //        printf("Hello World!\n");
    
63. 
    
64.         cout << "Call " << HugeCalc::GetVer() << endl << endl;
    
65. 
    
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));
    
80. 
    
81.     initsp();
    
82. 
    
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
    
93. 
    
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.         }
     
111.         
     
112.     }
     
113. 
     
114.         CONST UINT32 u32Timer_Filter = HugeCalc::GetTimer( TIMER_UNIT_ms );
     
115.     cout << "Filter cost "<<u32Timer_Filter<<"ms"<<endl;
     
116.         HugeCalc::ResetTimer();
     
117.     hugeStart-=4;
     
118. 
     
119.     int tcount=0,gcount=0;
     
120. 
     
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
     
142. 
     
143.         }
     
144.         hugeStart+=30;
     
145.     }
     
146. 
     
147. 
     
148.     // 记录计算耗时
     
149.         CONST UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_ms );
     
150. 
     
151. 
     
152.         // 输出
     
153.         cout << "Total took " << u32Timer_Calc << " ms" << endl;
     
154.     cout << "Total test "<<tcount<<" Numbers, get "<<gcount<<" numbers"<<endl;
     
155. 
     
156.         system( "pause" );
     
157. 
     
158.         return 0;
     
159. }

复制代码

gxqcn

原帖由 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

修改BUG后的程序运行数次分别得到数目为
15,0,16,18,16
所以数目同估计公式还是基本匹配的

mathe

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

我将19#重新编辑了下。

将 mathe 代码中的 initsp() 替换成 "HugeCalc::GetPrimeList( sptable, spcount, 0, SMALLPRIMES );"
并将所有 int 型变量替换成无符号型；以及其它一些小改动。。。

mathe

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

因为我们是从 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%