素数分解

无心人 · 发表于 2008-6-27 13:10:33

呵呵那复杂度也是迅速增加的我不知道到底多小合适

mathe · 发表于 2008-6-28 06:53:16

昨天重新从和数$1.8*10^15$左右开始搜索，得到一个满足除了条件a)的解，也就是同时满足本身素数，17个素数和，45个素数和，979个素数和以及2007个素数和的条件的: 1814523290364469

平均来说，应该每50个左右这样的数可以找到一个3个素数和的情况。如果找一个花费1天的时间，那么需要50天（或者100台机器花费半天）

无心人 · 发表于 2008-6-28 08:23:46

并行应该没这么大加速比吧

mathe · 发表于 2008-6-28 09:12:19

还好，这个问题并行性很好，而且复杂度不算太高。100台机器运行半天只是有可能产生结果，但是不一定产生，如果运行时间更加长一些，比如一星期，出结果的可能性就非常大了。比如假设运行一天没找到结果概率为$e^-1$,那么运行一星期没有找到结果的概率就是$e^-7=0.1%$

mathe · 发表于 2008-6-28 09:23:23

刚刚查看一下,又产生了一个 1852547403117883 下面是完整的代码,其中START_V可以修改成任何小于$2*10^18$,从而程序会从不同位置开始运行

// 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" )
#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 ((UINT64)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 979
#define MS 45
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;
}
#define START_V 1800000000000000LL
int main(int argc, char* argv[])
{
int i;
long long sumMF=0;
long long sumMS=0;
long long max_sum=2000000000000000000LL;
initSift(10000000000000000LL);
pMSblock = START_V/MS;
pMFblock = START_V/MF;
pMFcount=generalSift(pMFblock,pMFblock+BLOCK_SIZE,primeTabMF,TAB_LEN);
pMScount=generalSift(pMSblock,pMSblock+BLOCK_SIZE,primeTabMS,TAB_LEN);
printf("Init finished\n");
for(i=0;i<MS;i++){
sumMS+=primeTabMS[i];
}
for(i=0;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;
}

复制代码

mathe · 发表于 2008-6-28 09:25:39

此外,素数筛选算法在这里非常重要,但是我挑选的不是liangbch提供的最快的版本,而是修改比较方便的版本,如果能够将它替换成更快的版本,应该还有更加好的效率

无心人 · 发表于 2008-6-28 09:32:56

这个筛能否扩展到超过64位？

无心人 · 发表于 2008-6-28 09:42:23

粗略的看了你的代码超过64位后由于分段小会遇到很多素数，其倍数根本是不在筛里要白白计算其倍数现在就有这个问题但浪费的时间不会很大但随之搜索范围的扩大，将越来越浪费

mathe · 发表于 2008-6-28 09:45:00

依赖于前面我们预先产生的素数表格的范围。在这个程序中，我预先产生的素数表格为$10^8$,也就是用筛选法，最多能够筛到$10^16$，而程序中长度为979和2007个连续素数和都用到这个筛选结果，所以这个和不能超过不超过$10^16$的最大的979个素数和，这应该是个接近$10^19$的范围，但是没有达到$10^19$ 但是在这个程序里面其实$4*10^18$作为上界应该还是安全的，但是更加大一些可能会不安全，因为我们必须保证两个这么大的数相加不会溢出

无心人 · 发表于 2008-6-28 09:52:04

Start_V改成输入每隔10秒输出当前块的首数字你程序搜索到$4 xx 10^18$ 要多长时间？

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

[擂台] 素数分解

评分