数学研发论坛

 找回密码
 欢迎注册
楼主: 无心人

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

[复制链接]
发表于 2008-2-14 16:57:50 | 显示全部楼层
奇怪,实际计算表明用TestPrimality(5)去预选,又慢了很多,郭老大能猜测到底是什么原因吗?感觉挺奇怪的
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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个数,结果就比较相近了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-2-14 18:23:26 | 显示全部楼层
发现我的代码的一处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
            }
        }
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-2-14 18:23:58 | 显示全部楼层
新的代码:
  1. // HugeCalcDemo.cpp : Defines the entry point for the console application.
  2. //

  3. //        Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
  4. //                Win32 Debug:        Debug Multithreaded DLL
  5. //                Win32 Release:        Multithreaded DLL

  6. #include <time.h>
  7. #include <iostream>
  8. using namespace std;
  9. #include "../../../../HugeCalc_API/CppAPI/Include/HugeCalc.h"        // 公共接口
  10. #include "../../../../HugeCalc_API/CppAPI/Include/HugeInt.h"        // 10进制系统
  11. #include "../../../../HugeCalc_API/CppAPI/Include/HugeIntX.h"        // 16进制系统

  12. #pragma message( "automatic link to ../../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
  13. #pragma comment( lib, "../../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )

  14. //        http://www.swox.com/gmp/#TRY
  15. //        1 + gcd(87324,78263148,7896) * (10^1989879887 mod 471!)

  16. #if 1
  17.         #define _Test_HI
  18.         #define integer CHugeInt
  19. #else
  20.         #define _Test_HX
  21.         #define integer CHugeIntX
  22. #endif

  23. #define SQRTSMALLPRIMES 4000
  24. #define SMALLPRIMES (SQRTSMALLPRIMES*SQRTSMALLPRIMES)
  25. int spcount;
  26. int sptable[SMALLPRIMES];

  27. #define LIMIT_SIZE 10000000
  28. #define RANGE (LIMIT_SIZE<<5)
  29. int mask[LIMIT_SIZE];
  30. #define SET_MASK(x)  (mask[(x)>>5]|=1<<((x)&31))
  31. #define IS_SET(x)    (mask[(x)>>5]&(1<<((x)&31)))

  32. void initsp()
  33. {
  34.     int i,j;
  35.     memset(sptable,-1,sizeof(sptable));
  36.     sptable[0]=sptable[1]=0;
  37.     for(i=2;i<SQRTSMALLPRIMES;i++){
  38.         if(sptable[i]){
  39.             for(j=i*i;j<SMALLPRIMES;j+=i){
  40.                 sptable[j]=0;
  41.             }
  42.         }
  43.     }
  44.     for(i=0,spcount=0;i<SMALLPRIMES;i++){
  45.         if(sptable[i]){
  46.             sptable[spcount++]=i;
  47.         }
  48.     }
  49. }


  50. int main(int argc, char* argv[])
  51. {
  52. //        printf("Hello World!\n");

  53.         cout << "Call " << HugeCalc::GetVer() << endl << endl;

  54.         if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel())
  55.         {
  56.                 cout << endl << "警告:您未通过 HugeCalc.dll 的许可认证!" \
  57.                         << endl << endl << "解决方案可选下列方案之一:" \
  58.                         << endl <<        "    一、请将本程序移动到“/CopyrightByGuoXianqiang/[../]”目录下运行;" \
  59.                         << endl <<        "    二、或请在 HugeCalc.chm 中进行注册(一劳永逸)。" \
  60.                         << endl << endl;
  61.                 system( "pause" );
  62.                 return (-1);
  63.         }
  64.         // 初始化
  65.         HugeCalc::EnableTimer( TRUE );
  66.         HugeCalc::ResetTimer();
  67.     HugeCalc::SeedRandom(time(NULL));

  68.     initsp();

  69.     integer hugeStart, hugeTemp;
  70. #ifdef _Test_HI
  71.     hugeStart.Random(100);
  72. #else
  73.     hugeStart.Random(335);
  74. #endif
  75.     int m=hugeStart%30;
  76.     m=15-m;
  77.     if(m<0)m+=15;
  78.     hugeStart+=m;///After the normalization, hugeStart%30==15

  79.     int filter_list[]={-4,-2,2,4};
  80.     int i;
  81.     for(i=3;i<spcount;i++){//skip prime 2,3,5
  82.         int p=sptable[i];
  83.         m=hugeStart%p;
  84.         int v=HugeCalc::InvertMod(30,p);
  85.         int k;
  86.         for(k=0;k<sizeof(filter_list)/sizeof(filter_list[0]);k++){
  87.             long long s=-m+filter_list[k];
  88.             while(s<0)s+=p;
  89.             s*=v;
  90.             s%=p;//hugeStart+s*30==filter_list[k] (mod p)
  91.             int h;
  92.             for(h=(int)s;h<RANGE;h+=p){
  93.                 SET_MASK(h);///Filter hugeStart+h*30
  94.             }
  95.         }
  96.         
  97.     }

  98.         CONST UINT32 u32Timer_Filter = HugeCalc::GetTimer( TIMER_UNIT_ms );
  99.     cout << "Filter cost "<<u32Timer_Filter<<"ms"<<endl;
  100.         HugeCalc::ResetTimer();
  101.     hugeStart-=4;

  102.     int tcount=0,gcount=0;

  103.     for(i=0;i<RANGE;i++){
  104.         if(!IS_SET(i)){///hugeStart+i*30 may be a candidate
  105.             tcount++;
  106. #if 0
  107.             if(((hugeTemp=hugeStart)).TestPrimality(5)!=COMPOSITE&&
  108.                 (hugeTemp+=2).TestPrimality(5)!=COMPOSITE&&
  109.                 (hugeTemp+=4).TestPrimality(5)!=COMPOSITE&&
  110.                 (hugeTemp+=2).TestPrimality(5)!=COMPOSITE){
  111. #endif
  112.                 if(((hugeTemp=hugeStart)).IsPrime()&&
  113.                     (hugeTemp+=2).IsPrime()&&
  114.                     (hugeTemp+=4).IsPrime()&&
  115.                     (hugeTemp+=2).IsPrime()){
  116.                         ///Output the candidate;
  117.                            CONST LPCTSTR lpShow = hugeStart.GetStr( FS_DEFAULT );
  118.                     gcount++;
  119.                     cout<<"Found n,n+2,n+6,n+8 where n is "<<lpShow <<endl;
  120.                 }
  121. #if 0
  122.             }
  123. #endif

  124.         }
  125.         hugeStart+=30;
  126.     }


  127.     // 记录计算耗时
  128.         CONST UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_ms );


  129.         // 输出
  130.         cout << "Total took " << u32Timer_Calc << " ms" << endl;
  131.     cout << "Total test "<<tcount<<" Numbers, get "<<gcount<<" numbers"<<endl;

  132.         system( "pause" );

  133.         return 0;
  134. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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) 是个不错的选择?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-2-14 18:35:52 | 显示全部楼层
修改BUG后的程序运行数次分别得到数目为
15,0,16,18,16
所以数目同估计公式还是基本匹配的
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-2-15 10:37:45 | 显示全部楼层
我将19#重新编辑了下。

将 mathe 代码中的 initsp() 替换成 "HugeCalc::GetPrimeList( sptable, spcount, 0, SMALLPRIMES );"
并将所有 int 型变量替换成无符号型;以及其它一些小改动。。。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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.   }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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%
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2019-8-26 00:22 , Processed in 0.048204 second(s), 15 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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