数学研发论坛

 找回密码
 欢迎注册
楼主: mathe

[擂台] 素数分解

[复制链接]
发表于 2008-6-26 09:02:12 | 显示全部楼层
原帖由 mathe 于 2008-6-26 08:32 发表
我的想法是实际上上面代码中只有俩部分计算需要素数列表,也就是计算sumMF和sumMS部分.与其开始保存了所有素数的列表,不如俩部分都动态生成,比如保存两个长度为1,000,000左右大小范围内的素数,分别在当前iMF和iMS指向 ...


也是的。
而且你现在用的是全局变量,我一看是大数组,就没敢去调试了。

由于本题与进制无关(“回文素数”就与进制相关),
所以其效率应该是比较好的。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-6-26 09:06:29 | 显示全部楼层
原帖由 mathe 于 2008-6-26 08:58 发表

说明一下:
i)开始我采用start.NextPrime(start);这种用法,不过有点不放心这种用法是否可以(毕竟输入输出在同一个地址),帮助文件里面没有说明.此外,帮助中应该还要说明函数NextPrime()会修改*this,然后返回*this;
...


在帮助文档 HugeCalc.chm::快速生成最邻近素数 Remarks 第4条中有交待:
这里的入参均声明成“CONST”,但与HugeCalc其它函数不同的是:“A.NextPrime( B )”的意义是生成比“B”更大的素数“A”,对入参的CONST限定仅表明对入参进行保护,程序内部特意没有做“A非B”的检查,以便用户可以通过反复调用“A.NextPrime( A )”而得到一个类似超大素数发生器的工具;对于函数“PreviousPrime(..)”亦类似


A.NextPrime( B ) 中,当“A非B”时,B保持为定值;否则,B(即A)才可能被修改。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-6-26 09:06:51 | 显示全部楼层
原帖由 gxqcn 于 2008-6-26 09:02 发表


也是的。
而且你现在用的是全局变量,我一看是大数组,就没敢去调试了。

由于本题与进制无关(“回文素数”就与进制相关),
所以其效率应该是比较好的。

其实只要有足够大虚拟内存就没有问题,这个大数组也就顺序访问两次
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-6-26 10:44:20 | 显示全部楼层
现在改写liangbch的一个分段筛选代码,可以开始搜索了:

  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. // 以下代码源自 mathe 的,略有修改

  7. #include <time.h>
  8. #include <iostream>
  9. #include <stdio.h>
  10. #include <list>
  11. using namespace std;
  12. #include "../../../HugeCalc_API/CppAPI/Include/HugeCalc.h"  // 公共接口
  13. #include "../../../HugeCalc_API/CppAPI/Include/HugeInt.h"   // 10进制系统
  14. #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" )

  17. #define SRC_NAME "flterm.txt"
  18. #if 1
  19.     #define _Test_HI
  20.     #define integer CHugeInt
  21. #else
  22.     #define _Test_HX
  23.     #define integer CHugeIntX
  24. #endif

  25. #define MAX_N  203280220
  26. unsigned int plist[MAX_N];

  27. #include "math.h"
  28. #include "memory.h"
  29. #include "windows.h"

  30. //**********************************************
  31. // n 不能大于L2 cache,否则效能很低
  32. DWORD Sift32(DWORD n,DWORD *primeTab)
  33. // 对0到n 之间(不包括n)的数进行筛选,返回质数的个数
  34. //**********************************************
  35. {
  36.         BYTE *pBuff= new BYTE[n];
  37.         if (pBuff==NULL)
  38.                 return 0;

  39.         DWORD i,sqrt_n,p;
  40.         DWORD pn;
  41.         //----------------------------
  42.         sqrt_n=(DWORD)(sqrt((double)n)+1);
  43.        
  44.         while (sqrt_n * sqrt_n >=n)
  45.                 sqrt_n--;
  46.        
  47.         memset(pBuff,1,n*sizeof(BYTE));
  48.         *pBuff=0;
  49.         *(pBuff+1)=0; //1不是质数
  50.         *(pBuff+2)=1; //2是质数

  51.         for ( i=4;i<n;i+=2) //筛去2的倍数,不包括2
  52.                 *(pBuff+i)=0;
  53.        
  54.         for (i=3;i<=sqrt_n;) //i 是质数
  55.         {
  56.                 for (p=i*i;p<n; p+=2*i)
  57.                         *(pBuff+p)=0; //筛去i 的奇数倍
  58.                 i+=2;
  59.                 while ( *(pBuff+i)==0 )
  60.                         i++; //前进到下一个质数
  61.         }
  62.         pn=0; //素数个数
  63.         for (i=2;i<n;i++)
  64.         {
  65.                 if (pBuff[i])
  66.                 {
  67.                         primeTab[pn++]=i;
  68.                 }
  69.         }
  70.         delete[] pBuff;
  71.         return pn;
  72. }

  73. #define L2_CACHE (2*1024*1024)
  74. static DWORD getSiftBuffLength()
  75. {
  76.         DWORD n=L2_CACHE * 3/4;
  77.         DWORD rem=n % 30;
  78.         if (rem !=0)
  79.                 n-=rem;
  80.         return n;   //n must times of 30
  81. }


  82. //**********************************************
  83. // copy prime from psrc to primeTab
  84. // mode can be 4 or 8
  85. // if mode==4, primeTab will be consider as DWORD array
  86. // if mode==8, primeTab will be conseder ad UINT64 array
  87. // return total prime Count
  88. //**********************************************
  89. int copyL1PrimeToArray(const DWORD *psrc,int srcCount,void *primeTab,int primeBuffLen, int mode)
  90. {
  91.         DWORD pr;
  92.         int i;
  93.         pr=0;
  94.         if (mode==4)
  95.         {
  96.                 DWORD *UINT32PrimeTab=(DWORD*)primeTab;
  97.                 for (i=0;i<srcCount && pr<primeBuffLen;i++)
  98.                 {
  99.                         UINT32PrimeTab[pr]=psrc[i];
  100.                         pr++;
  101.                 }
  102.                 for (;i<srcCount;i++ )
  103.                 {        pr++; }
  104.         }
  105.         else if (mode==8)
  106.         {
  107.                 UINT64 *UINT64PrimeTab=(UINT64*)primeTab;
  108.                 for (i=0;i<srcCount && pr<primeBuffLen;i++)
  109.                 {
  110.                         UINT64PrimeTab[pr]=(UINT64)psrc[i];
  111.                         pr++;
  112.                 }
  113.                 for (;i<srcCount;i++ )
  114.                 {        pr++; }
  115.         }
  116.         return pr;
  117. }
  118. //**********************************************
  119. // mode can be 4 or 8
  120. // if mode==4, primeTab will be consider as DWORD array
  121. // if mode==8, primeTab will be conseder ad UINT64 array
  122. // return total prime Count
  123. //**********************************************
  124. static UINT64 GetPrimeFromSiftBuff(const BYTE *siftBuff,UINT64 begin, UINT64 end,
  125.                                            void *primeTab,UINT64 currCount,UINT64 primeBuffLen, int mode)
  126. {
  127.         UINT64 pr;
  128.         int i;
  129.         pr=currCount;
  130.         DWORD *UINT32PrimeTab=(DWORD*)primeTab;
  131.         UINT64 *UINT64PrimeTab=(UINT64*)primeTab;
  132.         //-----------------------------------
  133.         if (mode==4)
  134.         {
  135.                 for (i=0;i<end-begin && pr < primeBuffLen;i++)
  136.                 {
  137.                         if (siftBuff[i]==1)
  138.                         {
  139.                                 UINT64 p=begin+(UINT64)i;
  140.                                 UINT32PrimeTab[pr]=(DWORD)p;
  141.                                 pr++;
  142.                         }
  143.                 }
  144.         }
  145.         else if (mode==8)
  146.         {
  147.                 for (i=0;i<end-begin && pr < primeBuffLen;i++)
  148.                 {
  149.                         if (siftBuff[i]==1)
  150.                         {
  151.                                 UINT64 p=begin+(UINT64)i;
  152.                                 UINT64PrimeTab[pr]=p;
  153.                                 pr++;
  154.                         }
  155.                 }
  156.         }
  157.        
  158.         for (;i<end-begin;i++)
  159.         {
  160.                 if (siftBuff[i]==1)
  161.                 {
  162.                         pr++;
  163.                 }
  164.         }
  165.         return pr;
  166. }

  167. DWORD L1PrimeBuffLen;
  168. DWORD *L1PrimeTab;
  169. DWORD L1PrimeCount;
  170. void initSift(UINT64 n)
  171. {
  172.     DWORD sqrt_n;
  173.     double f;
  174.     sqrt_n=(DWORD)(sqrt((double)n));
  175.     sqrt_n=(sqrt_n+29)/30*30;
  176.     f=(double)sqrt_n*10.0/9.0+32;
  177.         L1PrimeBuffLen= (DWORD)f;
  178.        
  179.         L1PrimeTab= new DWORD[L1PrimeBuffLen];

  180.         L1PrimeCount= Sift32(sqrt_n,L1PrimeTab);
  181. }

  182. //**********************************************
  183. // 使用普通分段筛的方法筛素数,
  184. // // mode can be 4 or 8
  185. // if mode==4, primeTab will be consider as DWORD array
  186. // if mode==8, primeTab will be conseder ad UINT64 array
  187. UINT64 generalSift(UINT64 begin, UINT64 end,void *primeTab, UINT64 primeTabLimtCount)
  188. // 对0到n 之间(不包括n)的数进行筛选,返回质数的个数
  189. //**********************************************
  190. {
  191.         DWORD  siftBuffLen;
  192.         UINT64  pr=0;

  193.         siftBuffLen=end-begin;
  194.         BYTE *siftBuff=new BYTE[siftBuffLen];
  195.         if (siftBuff==NULL)
  196.                 return 0;
  197.        
  198.         {
  199.                 DWORD dwPrime,endPrime;
  200.                 UINT64 begIndex;
  201.                 DWORD  endIdx;
  202.                 int idx,i;
  203.                
  204.                 endIdx=end-begin;

  205.                 memset(siftBuff,1,siftBuffLen);
  206.                 endPrime= (DWORD)sqrt((double)end)+1;
  207.                
  208.                 while (endPrime*endPrime>=end)
  209.                         endPrime--;
  210.                 i=0;
  211.                 while (true)
  212.                 {
  213.                         DWORD inc;
  214.                         UINT64 t;
  215.                         dwPrime=L1PrimeTab[i];
  216.                         if (dwPrime * dwPrime>=end)
  217.                                 break;
  218.                         t=(begin+dwPrime-1)/dwPrime;
  219.                         if ( t<dwPrime)
  220.                                 t=dwPrime;
  221.                        
  222.                         if (dwPrime==2)
  223.                         {
  224.                                 idx=(int)(t * dwPrime-begin);
  225.                                 while (idx<endIdx)
  226.                                 {
  227.                                         siftBuff[idx]=0;
  228.                                         idx+=dwPrime;
  229.                                 }
  230.                         }
  231.                         else
  232.                         {
  233.                                 if ( t % 2==0)
  234.                                         t++;                //素数的奇数倍
  235.                                 inc= dwPrime+dwPrime;
  236.                                 idx=(int)(t * dwPrime-begin);
  237.                                 while (idx<endIdx)
  238.                                 {
  239.                                         siftBuff[idx]=0;
  240.                                         idx+=inc;
  241.                                 }
  242.                         }
  243.                         i++;
  244.                 }
  245.                 pr=GetPrimeFromSiftBuff(siftBuff,begin,end,primeTab,pr,primeTabLimtCount,8);
  246.         }
  247.         delete[] siftBuff;
  248.         return pr;
  249. }


  250. int search_match_prime(unsigned value)
  251. {
  252.     int start=0,mid,end=MAX_N-1;
  253.     mid=(start+end)/2;
  254.     while(end-start>=2){
  255.         if(plist[mid]<value){
  256.             start=mid;
  257.         }else{
  258.             end=mid;
  259.         }
  260.         mid=(start+end)/2;
  261.     }
  262.     return mid;
  263. }

  264. bool is_list(int num, unsigned long long sum)
  265. {
  266.     unsigned long long aver=sum/num;
  267.     integer start(aver),end,sump(0);
  268.     start.NextPrime(start);
  269.     int prevc,nextc;
  270.     prevc=num/2;
  271.     nextc=num-1-prevc;
  272.     int i;
  273.     sump=end=start;
  274.     for(i=0;i<nextc;i++){
  275.         integer tmp(end);
  276.         end.NextPrime(tmp);
  277.         sump+=end;
  278.     }
  279.     for(i=0;i<prevc;i++){
  280.         integer tmp(start);
  281.         start.PreviousPrime(tmp);
  282.         sump+=start;
  283.     }
  284.     if(sump<integer(sum)){
  285.         while(sump<integer(sum)){
  286.             integer newp;
  287.             newp.NextPrime(end);
  288.             sump+=newp-start;
  289.             integer tmp(start);
  290.             start.NextPrime(tmp);
  291.             end=newp;
  292.         }
  293.     }else{
  294.         while(sump>integer(sum)){
  295.             integer newp;
  296.             newp.PreviousPrime(start);
  297.             sump-=end-newp;
  298.             integer tmp(end);
  299.             end.PreviousPrime(tmp);
  300.             start=newp;
  301.         }
  302.     }
  303.     return sump==integer(sum);
  304. }

  305. #define MF  2007
  306. #define MS  979

  307. void candidate(unsigned long long x)
  308. {
  309.     integer y((unsigned long long)x);
  310.     if(!y.IsPrime())
  311.         return;
  312.     if(!is_list(45,x))
  313.         return;
  314.     printf("Cand45: %lld\n",x);
  315.     if(!is_list(17,x))
  316.         return;
  317.     printf("Cand17: %lld\n",x);
  318.     if(!is_list(3,x))
  319.         return;
  320.     printf("Find: %lld\n",x);
  321. }

  322. #define BLOCK_SIZE 10000000
  323. #define TAB_LEN 689388
  324. UINT64 primeTabMF[MF+TAB_LEN];
  325. UINT64 primeTabMS[MS+TAB_LEN];
  326. UINT64 pMFblock,pMSblock;
  327. UINT64 pMFcount,pMScount;
  328. void msNextBlock()
  329. {
  330.     memmove(primeTabMS,primeTabMS+pMScount-MS+1,(MS-1)*sizeof(UINT64));
  331. start:
  332.     pMSblock+=BLOCK_SIZE;
  333.     pMScount=generalSift(pMSblock,pMSblock+BLOCK_SIZE,primeTabMS+MS-1,TAB_LEN);
  334.     if(pMScount==0)goto start;
  335.     pMScount+=MS-1;
  336. }

  337. void mfNextBlock()
  338. {
  339.     memmove(primeTabMF,primeTabMF+pMFcount-MF+1,(MF-1)*sizeof(UINT64));   
  340. start:
  341.     pMFblock+=BLOCK_SIZE;
  342.     pMFcount=generalSift(pMFblock,pMFblock+BLOCK_SIZE,primeTabMF+MF-1,TAB_LEN);
  343.     if(pMFcount==0)goto start;
  344.     pMFcount+=MF-1;
  345. }

  346. int main(int argc, char* argv[])
  347. {
  348.     int i;
  349.     long long sumMF=0;
  350.     long long sumMS=0;
  351.     long long max_sum=2000000000000000000LL;
  352.     initSift(10000000000000000LL);
  353.     pMFcount=generalSift(0,BLOCK_SIZE,primeTabMF,TAB_LEN);
  354.     pMScount=generalSift(0,BLOCK_SIZE,primeTabMS,TAB_LEN);

  355.     printf("Init finished\n");
  356.     for(i=0;i<MS;i++){
  357.         sumMS+=primeTabMF[i];
  358.     }
  359.     sumMF=sumMS;
  360.     for(;i<MF;i++){
  361.         sumMF+=primeTabMF[i];
  362.     }
  363.     int iMF=0,iMS=0;
  364.     while(sumMF<=max_sum){
  365.         while(sumMS<sumMF){
  366.             sumMS+=primeTabMS[iMS+MS]-primeTabMS[iMS];
  367.             iMS++;
  368.             if(iMS>=pMScount-MS){
  369.                 msNextBlock();
  370.                 iMS=0;
  371.             }
  372.         }
  373.         if(sumMS==sumMF){
  374.             candidate(sumMF);
  375.         }
  376.         sumMF+=primeTabMF[iMF+MF]-primeTabMF[iMF];
  377.         iMF++;
  378.         if(iMF>=pMFcount-MF){
  379.             mfNextBlock();
  380.             iMF=0;
  381.         }
  382.     }
  383.     return 0;
  384. }
复制代码
不知道代码是否正确。
现在只找到两个满足本身素数以及表示成45,979,2007个素数和的结果:

Cand45: 1453832430671
Cand45: 3937966662683
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-6-26 10:53:58 | 显示全部楼层
程序还有bug,运行到某一步会出错.不过这个程序调试起来太麻烦了,要运行一段时间才能够出错
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-6-26 11:14:34 | 显示全部楼层
找到一个bug,
generalshft函数中:
                        if (dwPrime * dwPrime>=end)
                                break;
改为:
                        if ((UINT64)dwPrime * dwPrime>=end)
                                break;
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-6-26 11:28:15 | 显示全部楼层
谁来运行一下看看?这个程序运行完了估计能够找出一个结果来。但是估计要运行数天。我的机器明天一大早要断电一次,肯定运行不完
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-6-26 11:51:22 | 显示全部楼层
你做好了,我帮你运行
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-6-26 13:55:02 | 显示全部楼层
就上面的代码,如36楼修改了一行就可以了。但是要调用gxqcn的HugeCalc.
我的机器现在计算到和为$2*10^15$左右,还是没有一个同时是2007,979,45,17个素数和的素数。现在打印出来都是到45个素数和为止,最后几个为:
Cand45: 183071317698557
Cand45: 201340432884221
Cand45: 206978108198867
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-6-26 15:03:03 | 显示全部楼层


你给我个编译后的程序,我服务器新做了
没装编译器

linux下倒可以尝试,呵呵
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2019-1-16 12:40 , Processed in 0.051645 second(s), 14 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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