奇怪，实际计算表明用TestPrimality(5)去预选，又慢了很多，郭老大能猜测到底是什么原因吗？感觉挺奇怪的
还有我上面测试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个数，结果就比较相近了。
发现我的代码的一处BUG.
下面代码中int s;应该改成long long s;
要不然在SQRTSMALLPRIMES>215时计算就有可能发生溢出，从而出现错误的计算，导致结果偏少。
修改以后，明显找到的结果要增加很多。
还有一点发现的结果是SQRTSMALLPRIMES的增加对速度影响不是很大，所以可以适当提高，我的机器上增加到4000后计算速度
可以有显著提高，总运算时间基本都在2分钟之内（筛选1分钟，后面测试1分钟）
for(k=0;k<sizeof(filter_list)/sizeof(filter_list);k++){
int s=-m+filter_list;
while(s<0)s+=p;
s*=v;
s%=p;//hugeStart+s*30==filter_list (mod p)
int h;
for(h=s;h<RANGE;h+=p){
SET_MASK(h);///Filter hugeStart+h*30
}
}
新的代码:// HugeCalcDemo.cpp : Defines the entry point for the console application.
//
// Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
// Win32 Debug: Debug Multithreaded DLL
// Win32 Release: Multithreaded DLL
#include <time.h>
#include <iostream>
using namespace std;
#include "../../../../HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
#include "../../../../HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
#include "../../../../HugeCalc_API/CppAPI/Include/HugeIntX.h" // 16进制系统
#pragma message( "automatic link to ../../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#pragma comment( lib, "../../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
// http://www.swox.com/gmp/#TRY
// 1 + gcd(87324,78263148,7896) * (10^1989879887 mod 471!)
#if 1
#define _Test_HI
#define integer CHugeInt
#else
#define _Test_HX
#define integer CHugeIntX
#endif
#define SQRTSMALLPRIMES 4000
#define SMALLPRIMES (SQRTSMALLPRIMES*SQRTSMALLPRIMES)
int spcount;
int sptable;
#define LIMIT_SIZE 10000000
#define RANGE (LIMIT_SIZE<<5)
int mask;
#define SET_MASK(x)(mask[(x)>>5]|=1<<((x)&31))
#define IS_SET(x) (mask[(x)>>5]&(1<<((x)&31)))
void initsp()
{
int i,j;
memset(sptable,-1,sizeof(sptable));
sptable=sptable=0;
for(i=2;i<SQRTSMALLPRIMES;i++){
if(sptable){
for(j=i*i;j<SMALLPRIMES;j+=i){
sptable=0;
}
}
}
for(i=0,spcount=0;i<SMALLPRIMES;i++){
if(sptable){
sptable=i;
}
}
}
int main(int argc, char* argv[])
{
// printf("Hello World!\n");
cout << "Call " << HugeCalc::GetVer() << endl << endl;
if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel())
{
cout << endl << "警告：您未通过 HugeCalc.dll 的许可认证！" \
<< endl << endl << "解决方案可选下列方案之一：" \
<< endl << " 一、请将本程序移动到“/CopyrightByGuoXianqiang/[../]”目录下运行；" \
<< endl << " 二、或请在 HugeCalc.chm 中进行注册（一劳永逸）。" \
<< endl << endl;
system( "pause" );
return (-1);
}
// 初始化
HugeCalc::EnableTimer( TRUE );
HugeCalc::ResetTimer();
HugeCalc::SeedRandom(time(NULL));
initsp();
integer hugeStart, hugeTemp;
#ifdef _Test_HI
hugeStart.Random(100);
#else
hugeStart.Random(335);
#endif
int m=hugeStart%30;
m=15-m;
if(m<0)m+=15;
hugeStart+=m;///After the normalization, hugeStart%30==15
int filter_list[]={-4,-2,2,4};
int i;
for(i=3;i<spcount;i++){//skip prime 2,3,5
int p=sptable;
m=hugeStart%p;
int v=HugeCalc::InvertMod(30,p);
int k;
for(k=0;k<sizeof(filter_list)/sizeof(filter_list);k++){
long long s=-m+filter_list;
while(s<0)s+=p;
s*=v;
s%=p;//hugeStart+s*30==filter_list (mod p)
int h;
for(h=(int)s;h<RANGE;h+=p){
SET_MASK(h);///Filter hugeStart+h*30
}
}
}
CONST UINT32 u32Timer_Filter = HugeCalc::GetTimer( TIMER_UNIT_ms );
cout << "Filter cost "<<u32Timer_Filter<<"ms"<<endl;
HugeCalc::ResetTimer();
hugeStart-=4;
int tcount=0,gcount=0;
for(i=0;i<RANGE;i++){
if(!IS_SET(i)){///hugeStart+i*30 may be a candidate
tcount++;
#if 0
if(((hugeTemp=hugeStart)).TestPrimality(5)!=COMPOSITE&&
(hugeTemp+=2).TestPrimality(5)!=COMPOSITE&&
(hugeTemp+=4).TestPrimality(5)!=COMPOSITE&&
(hugeTemp+=2).TestPrimality(5)!=COMPOSITE){
#endif
if(((hugeTemp=hugeStart)).IsPrime()&&
(hugeTemp+=2).IsPrime()&&
(hugeTemp+=4).IsPrime()&&
(hugeTemp+=2).IsPrime()){
///Output the candidate;
CONST LPCTSTR lpShow = hugeStart.GetStr( FS_DEFAULT );
gcount++;
cout<<"Found n,n+2,n+6,n+8 where n is "<<lpShow <<endl;
}
#if 0
}
#endif
}
hugeStart+=30;
}
// 记录计算耗时
CONST UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_ms );
// 输出
cout << "Total took " << u32Timer_Calc << " ms" << endl;
cout << "Total test "<<tcount<<" Numbers, get "<<gcount<<" numbers"<<endl;
system( "pause" );
return 0;
}
原帖由 mathe 于 2008-2-14 16:57 发表 http://bbs.emath.ac.cn/images/common/back.gif
奇怪，实际计算表明用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) 是个不错的选择？
修改BUG后的程序运行数次分别得到数目为
15,0,16,18,16
所以数目同估计公式还是基本匹配的
Faint:if(m<0)m+=15;shoulbe be changed to=>if(m<0)m+=30;That's why there're so much zeros.
我将19#重新编辑了下。
将 mathe 代码中的 initsp() 替换成 "HugeCalc::GetPrimeList( sptable, spcount, 0, SMALLPRIMES );"
并将所有 int 型变量替换成无符号型；以及其它一些小改动。。。
Another small optimization opportunity:
change loopfor(i=0;i<RANGE;i++){
if(!IS_SET(i)){//hugeStart+i*30 may be a candidate
...
}
}tofor(i=0;i<LIMIT_SIZE;i++){
if(!mask){
hugeStart+=30*32;
}else{
for(t=0;t<32;t++){
if(mask&(1<<t)){
.....
}
hugeStart+=30;
}
}
}
因为我们是从 p --> p+2 --> p+6 --> p+8 这个顺序逐步测试素性的，应从“粗放”到“精细”过渡，
所以排在前面的数应以淘汰效率为主；排在后面的数应以准确度为主，比方说最后一个数就应直接调用一次 IsPrime() 即可。
我将原来的代码：if( (hugeTest=hugeStart).IsPrime()&&
(hugeTest+=2).IsPrime()&&
(hugeTest+=4).IsPrime()&&
(hugeTest+=2).IsPrime()){
//Output the candidate;替换成新的（已在19#楼修正）：if( (hugeTest=hugeStart).TestPrimality(1) != COMPOSITE &&
(hugeTest+=2).TestPrimality(2) != COMPOSITE &&
(hugeTest+=4).TestPrimality(3) != COMPOSITE &&
(hugeTest+=2).IsPrime() &&
(hugeTest-=2).IsPrime() &&
(hugeTest-=4).IsPrime() &&
(hugeTest-=2).IsPrime()){
//Output the candidate;搜索10^100后的连续13组“四生素数”，时间由先前的1'22.174''降为1’13.607''，提速约11.64%
补记：19:00 左右，再次对代码优化，去掉了临时变量“hugeTest”，节省了一次大数拷贝动作。但却用了一个C/C++少有的语法——在if逻辑语句中大量使用逗号“,”运算操作符。经测试，可在新标准下再提速18.31%