素数分解

gxqcn · 发表于 2008-6-26 09:02:12

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

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

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

gxqcn · 发表于 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)才可能被修改。

mathe · 发表于 2008-6-26 09:06:51

原帖由 gxqcn 于 2008-6-26 09:02 发表

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

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

其实只要有足够大虚拟内存就没有问题,这个大数组也就顺序访问两次

mathe · 发表于 2008-6-26 10:44:20

现在改写liangbch的一个分段筛选代码，可以开始搜索了:

// 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
// 以下代码源自 mathe 的，略有修改
#include <time.h>
#include <iostream>
#include <stdio.h>
#include <list>
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" )
#define SRC_NAME "flterm.txt"
#if 1
#define _Test_HI
#define integer CHugeInt
#else
#define _Test_HX
#define integer CHugeIntX
#endif
#define MAX_N 203280220
unsigned int plist[MAX_N];
#include "math.h"
#include "memory.h"
#include "windows.h"
//**********************************************
// n 不能大于L2 cache，否则效能很低
DWORD Sift32(DWORD n,DWORD *primeTab)
// 对0到n 之间(不包括n)的数进行筛选，返回质数的个数
//**********************************************
{
BYTE *pBuff= new BYTE[n];
if (pBuff==NULL)
return 0;
DWORD i,sqrt_n,p;
DWORD pn;
//----------------------------
sqrt_n=(DWORD)(sqrt((double)n)+1);
while (sqrt_n * sqrt_n >=n)
sqrt_n--;
memset(pBuff,1,n*sizeof(BYTE));
*pBuff=0;
*(pBuff+1)=0; //1不是质数
*(pBuff+2)=1; //2是质数
for ( i=4;i<n;i+=2) //筛去2的倍数，不包括2
*(pBuff+i)=0;
for (i=3;i<=sqrt_n;) //i 是质数
{
for (p=i*i;p<n; p+=2*i)
*(pBuff+p)=0; //筛去i 的奇数倍
i+=2;
while ( *(pBuff+i)==0 )
i++; //前进到下一个质数
}
pn=0; //素数个数
for (i=2;i<n;i++)
{
if (pBuff[i])
{
primeTab[pn++]=i;
}
}
delete[] pBuff;
return pn;
}
#define L2_CACHE (2*1024*1024)
static DWORD getSiftBuffLength()
{
DWORD n=L2_CACHE * 3/4;
DWORD rem=n % 30;
if (rem !=0)
n-=rem;
return n; //n must times of 30
}
//**********************************************
// copy prime from psrc to primeTab
// mode can be 4 or 8
// if mode==4, primeTab will be consider as DWORD array
// if mode==8, primeTab will be conseder ad UINT64 array
// return total prime Count
//**********************************************
int copyL1PrimeToArray(const DWORD *psrc,int srcCount,void *primeTab,int primeBuffLen, int mode)
{
DWORD pr;
int i;
pr=0;
if (mode==4)
{
DWORD *UINT32PrimeTab=(DWORD*)primeTab;
for (i=0;i<srcCount && pr<primeBuffLen;i++)
{
UINT32PrimeTab[pr]=psrc[i];
pr++;
}
for (;i<srcCount;i++ )
{ pr++; }
}
else if (mode==8)
{
UINT64 *UINT64PrimeTab=(UINT64*)primeTab;
for (i=0;i<srcCount && pr<primeBuffLen;i++)
{
UINT64PrimeTab[pr]=(UINT64)psrc[i];
pr++;
}
for (;i<srcCount;i++ )
{ pr++; }
}
return pr;
}
//**********************************************
// mode can be 4 or 8
// if mode==4, primeTab will be consider as DWORD array
// if mode==8, primeTab will be conseder ad UINT64 array
// return total prime Count
//**********************************************
static UINT64 GetPrimeFromSiftBuff(const BYTE *siftBuff,UINT64 begin, UINT64 end,
void *primeTab,UINT64 currCount,UINT64 primeBuffLen, int mode)
{
UINT64 pr;
int i;
pr=currCount;
DWORD *UINT32PrimeTab=(DWORD*)primeTab;
UINT64 *UINT64PrimeTab=(UINT64*)primeTab;
//-----------------------------------
if (mode==4)
{
for (i=0;i<end-begin && pr < primeBuffLen;i++)
{
if (siftBuff[i]==1)
{
UINT64 p=begin+(UINT64)i;
UINT32PrimeTab[pr]=(DWORD)p;
pr++;
}
}
}
else if (mode==8)
{
for (i=0;i<end-begin && pr < primeBuffLen;i++)
{
if (siftBuff[i]==1)
{
UINT64 p=begin+(UINT64)i;
UINT64PrimeTab[pr]=p;
pr++;
}
}
}
for (;i<end-begin;i++)
{
if (siftBuff[i]==1)
{
pr++;
}
}
return pr;
}
DWORD L1PrimeBuffLen;
DWORD *L1PrimeTab;
DWORD L1PrimeCount;
void initSift(UINT64 n)
{
DWORD sqrt_n;
double f;
sqrt_n=(DWORD)(sqrt((double)n));
sqrt_n=(sqrt_n+29)/30*30;
f=(double)sqrt_n*10.0/9.0+32;
L1PrimeBuffLen= (DWORD)f;
L1PrimeTab= new DWORD[L1PrimeBuffLen];
L1PrimeCount= Sift32(sqrt_n,L1PrimeTab);
}
//**********************************************
// 使用普通分段筛的方法筛素数,
// // mode can be 4 or 8
// if mode==4, primeTab will be consider as DWORD array
// if mode==8, primeTab will be conseder ad UINT64 array
UINT64 generalSift(UINT64 begin, UINT64 end,void *primeTab, UINT64 primeTabLimtCount)
// 对0到n 之间(不包括n)的数进行筛选，返回质数的个数
//**********************************************
{
DWORD siftBuffLen;
UINT64 pr=0;
siftBuffLen=end-begin;
BYTE *siftBuff=new BYTE[siftBuffLen];
if (siftBuff==NULL)
return 0;
{
DWORD dwPrime,endPrime;
UINT64 begIndex;
DWORD endIdx;
int idx,i;
endIdx=end-begin;
memset(siftBuff,1,siftBuffLen);
endPrime= (DWORD)sqrt((double)end)+1;
while (endPrime*endPrime>=end)
endPrime--;
i=0;
while (true)
{
DWORD inc;
UINT64 t;
dwPrime=L1PrimeTab[i];
if (dwPrime * dwPrime>=end)
break;
t=(begin+dwPrime-1)/dwPrime;
if ( t<dwPrime)
t=dwPrime;
if (dwPrime==2)
{
idx=(int)(t * dwPrime-begin);
while (idx<endIdx)
{
siftBuff[idx]=0;
idx+=dwPrime;
}
}
else
{
if ( t % 2==0)
t++; //素数的奇数倍
inc= dwPrime+dwPrime;
idx=(int)(t * dwPrime-begin);
while (idx<endIdx)
{
siftBuff[idx]=0;
idx+=inc;
}
}
i++;
}
pr=GetPrimeFromSiftBuff(siftBuff,begin,end,primeTab,pr,primeTabLimtCount,8);
}
delete[] siftBuff;
return pr;
}
int search_match_prime(unsigned value)
{
int start=0,mid,end=MAX_N-1;
mid=(start+end)/2;
while(end-start>=2){
if(plist[mid]<value){
start=mid;
}else{
end=mid;
}
mid=(start+end)/2;
}
return mid;
}
bool is_list(int num, unsigned long long sum)
{
unsigned long long aver=sum/num;
integer start(aver),end,sump(0);
start.NextPrime(start);
int prevc,nextc;
prevc=num/2;
nextc=num-1-prevc;
int i;
sump=end=start;
for(i=0;i<nextc;i++){
integer tmp(end);
end.NextPrime(tmp);
sump+=end;
}
for(i=0;i<prevc;i++){
integer tmp(start);
start.PreviousPrime(tmp);
sump+=start;
}
if(sump<integer(sum)){
while(sump<integer(sum)){
integer newp;
newp.NextPrime(end);
sump+=newp-start;
integer tmp(start);
start.NextPrime(tmp);
end=newp;
}
}else{
while(sump>integer(sum)){
integer newp;
newp.PreviousPrime(start);
sump-=end-newp;
integer tmp(end);
end.PreviousPrime(tmp);
start=newp;
}
}
return sump==integer(sum);
}
#define MF 2007
#define MS 979
void candidate(unsigned long long x)
{
integer y((unsigned long long)x);
if(!y.IsPrime())
return;
if(!is_list(45,x))
return;
printf("Cand45: %lld\n",x);
if(!is_list(17,x))
return;
printf("Cand17: %lld\n",x);
if(!is_list(3,x))
return;
printf("Find: %lld\n",x);
}
#define BLOCK_SIZE 10000000
#define TAB_LEN 689388
UINT64 primeTabMF[MF+TAB_LEN];
UINT64 primeTabMS[MS+TAB_LEN];
UINT64 pMFblock,pMSblock;
UINT64 pMFcount,pMScount;
void msNextBlock()
{
memmove(primeTabMS,primeTabMS+pMScount-MS+1,(MS-1)*sizeof(UINT64));
start:
pMSblock+=BLOCK_SIZE;
pMScount=generalSift(pMSblock,pMSblock+BLOCK_SIZE,primeTabMS+MS-1,TAB_LEN);
if(pMScount==0)goto start;
pMScount+=MS-1;
}
void mfNextBlock()
{
memmove(primeTabMF,primeTabMF+pMFcount-MF+1,(MF-1)*sizeof(UINT64));
start:
pMFblock+=BLOCK_SIZE;
pMFcount=generalSift(pMFblock,pMFblock+BLOCK_SIZE,primeTabMF+MF-1,TAB_LEN);
if(pMFcount==0)goto start;
pMFcount+=MF-1;
}
int main(int argc, char* argv[])
{
int i;
long long sumMF=0;
long long sumMS=0;
long long max_sum=2000000000000000000LL;
initSift(10000000000000000LL);
pMFcount=generalSift(0,BLOCK_SIZE,primeTabMF,TAB_LEN);
pMScount=generalSift(0,BLOCK_SIZE,primeTabMS,TAB_LEN);
printf("Init finished\n");
for(i=0;i<MS;i++){
sumMS+=primeTabMF[i];
}
sumMF=sumMS;
for(;i<MF;i++){
sumMF+=primeTabMF[i];
}
int iMF=0,iMS=0;
while(sumMF<=max_sum){
while(sumMS<sumMF){
sumMS+=primeTabMS[iMS+MS]-primeTabMS[iMS];
iMS++;
if(iMS>=pMScount-MS){
msNextBlock();
iMS=0;
}
}
if(sumMS==sumMF){
candidate(sumMF);
}
sumMF+=primeTabMF[iMF+MF]-primeTabMF[iMF];
iMF++;
if(iMF>=pMFcount-MF){
mfNextBlock();
iMF=0;
}
}
return 0;
}

复制代码

不知道代码是否正确。
现在只找到两个满足本身素数以及表示成45,979,2007个素数和的结果:

Cand45: 1453832430671
Cand45: 3937966662683

mathe · 发表于 2008-6-26 10:53:58

程序还有bug,运行到某一步会出错.不过这个程序调试起来太麻烦了,要运行一段时间才能够出错

mathe · 发表于 2008-6-26 11:14:34

找到一个bug,
generalshft函数中:
                     if (dwPrime * dwPrime>=end)
                              break;
改为:
                     if ((UINT64)dwPrime * dwPrime>=end)
                              break;

mathe · 发表于 2008-6-26 11:28:15

谁来运行一下看看？这个程序运行完了估计能够找出一个结果来。但是估计要运行数天。我的机器明天一大早要断电一次，肯定运行不完

无心人 · 发表于 2008-6-26 11:51:22

你做好了，我帮你运行

mathe · 发表于 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下倒可以尝试，呵呵

账号		自动登录	找回密码
密码			欢迎注册

[擂台] 素数分解