素数分解

无心人 · 发表于 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$
要多长时间？

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

[擂台] 素数分解

评分