K-Smith numbers

无心人 · 发表于 2008-11-25 18:57:55

ksmith100000000.TXT (46.03 KB, 下载次数: 3)

无心人 · 发表于 2008-11-25 18:59:10

1小时14分钟

liangbch · 发表于 2008-11-26 09:29:37

看来无心人你这个程序还是难敌C程序呀。到今日凌晨5点，我也完成了我的C程序。这个程序的复杂度近似于O(n)，在PM 1.7上,计算10^8以内的smith数(不输出)，不到10秒，好像是4.2秒，可是U盘忘带了。等明天吧，明天我贴出程序。

无心人 · 发表于 2008-11-26 09:40:44

不是不敌
是我初学
写不出高效的算法啊

如果C用的时间是1
Haskell通常在2-5

无心人 · 发表于 2008-11-26 15:00:59

用筛法做，是最佳的选择了，呵呵
不过啊
haskell实现筛法，似乎有点挑战啊
不是haskell不行
是俺笨哦

liangbch · 发表于 2008-11-26 15:42:14

使用晒法，我在5楼就提到了，关于复杂度我也在5楼提到了。只是在实现的时候发现不需使用内存池和链表。
另外，可以采用一些优化手段。我的做法是事先算出10000以下所有的数的数字和，用的时候查表。

无心人 · 发表于 2008-11-26 17:14:39

呵呵

这里我并不奢望多高等的算法
当学习haskell了
目前的估计是
30分钟
用的是逐个分解的算法
另外，算100000000个数字的数字和
是5秒

liangbch · 发表于 2008-11-27 10:01:45

贴出我的解法：
  1.欲求n以内的所有smith数，首先需要求出sqrt(n)以内的所有素数。
  2.为了快速的求出数x的数字和，事先计算出10-10000的数字和，在需要的时候直接查表。
  3.需要对 n以内的每个数分解质因数，对x分解质因数，最笨的方法是用sqrt(x)以内的素数逐个试除，
我这里采用类似用筛法求素数的方法，平均下来，对n以内的每个数分解质因数的复杂度为log(log(n).

#include "stdlib.h"
#include "stdio.h"
#include "math.h"
#include "memory.h"
typedef unsigned __int64 UINT64;
typedef unsigned long DWORD;
typedef unsigned short WORD;
typedef unsigned char BYTE;
//文件描述：使用最简单的方法筛素数
// 一定范围内质数的个数
// 10 4
// 100 25
// 1,000 168
// 10,000 1,229
// 100,000 9,592
// 1,000,000 78,498
// 10,000,000 664,579
// 100,000,000 5,761,455
// 1,000,000,000 50,847,534
// 10,000,000,000 455,052,511
// 100,000,000,000 4,118,054,813
// 1,000,000,000,000 37,607,912,018
// 10,000,000,000,000 346,065,536,839
// 100,000,000,000,000 3,204,941,750,802
// 1,000,000,000,000,000 29,844,570,422,669
// 10,000,000,000,000,000 279,238,341,033,925
// 100,000,000,000,000,000 2,623,557,157,654,233
// 1,000,000,000,000,000,000 24,739,954,287,740,860
// 估计x 范围内的质数的个数，欧拉公式 pi(x)=x /( ln(x));
// 经过测试，
// 在1-1百万之间，相邻2个质数的最大差为114,这两个质数分别为492113,492227
// 在1-1千万之间，相邻2个质数的最大差为154,这两个质数分别为4652353,4652507
// 在1-3千万之间，相邻2个质数的最大差为210,这两个质数分别为20831323,20831533
// 在1-1亿之间，相邻2个质数的最大差为220,这两个质数分别为47326693,47326913
//计算n以内素数个数的上限
UINT64 maxPrimeCount(UINT64 n)
{
double f= (double)n / log((double)n) * 10.0/9.0 +32;
return (UINT64)f;
}
DWORD Sift32(DWORD n,DWORD *primeTab)
// 对0到n 之间(不包括n)的数进行筛选，返回质数的个数
// n 不能大于L2 cache，否则效能很低
//**********************************************
{
BYTE *pBuff= new BYTE[n];
if (pBuff==NULL)
return 0;
DWORD i,sqrt_n,p;
DWORD pn;
//----------------------------
sqrt_n=(DWORD)(sqrt(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;
}

复制代码

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <memory.h>
#include <windows.h>
#include <assert.h>
typedef unsigned __int64 UINT64;
typedef unsigned long DWORD;
typedef unsigned short WORD;
typedef unsigned char BYTE;
#define SIFT_BLOCK_LEN (64*1024)
typedef struct fac_struct
{
DWORD n;
WORD digSum;
}FAC_ST;
BYTE g_digSumTab[10000];
FAC_ST g_facArray[SIFT_BLOCK_LEN+1];
extern UINT64 maxPrimeCount(UINT64 n);
extern DWORD Sift32(DWORD n,DWORD *primeTab);
WORD getDigSum(DWORD n)
{
DWORD s=0;
while (n>0)
{
s+= (n % 10);
n/=10;
}
return (WORD)s;
}
WORD getDigSum2(DWORD n)
{
assert(g_digSumTab[1]==1);
DWORD s=0;
if (n<10000)
return g_digSumTab[n];
else if (n<100000000)
{
DWORD t=n % 10000;
n/=10000;
return g_digSumTab[n]+g_digSumTab[t];
}
else
{
DWORD t0,t1;
t0=n % 10000;
n/=10000;
t1=n % 10000;
n/=10000;
return g_digSumTab[n]+g_digSumTab[t0]+g_digSumTab[t1];
}
}
void initTable(BYTE *tab,DWORD n)
{
DWORD i;
assert(n==10 || n==100 || n==1000 || n==10000);
for (i=0;i<n;i++)
{
tab[i]=(BYTE)getDigSum(i);
}
}
void searchSmithNumber(DWORD n,const char *fileName)
{
DWORD *primeTab=NULL;
DWORD base;
DWORD s_n;
DWORD primeCount;
DWORD i,count;
DWORD time1,time2;
DWORD c;
DWORD numBuff[1024];
int countInBuff;
FILE *fp=NULL;
if (fileName!=NULL)
{
fp=fopen(fileName,"wb");
}
time1=GetTickCount();
initTable(g_digSumTab,10000 );
s_n =(DWORD)sqrt((double)n)+2;
while (s_n * s_n >n)
s_n--;
//估算s_n以内素数个数的上限;
c=(DWORD)maxPrimeCount(s_n);
primeTab=new DWORD[c];
if (primeTab==NULL )
{
goto thisExit;
}
primeCount=Sift32(s_n+1,primeTab);
count=0; base=0;
countInBuff=0;
while (base<=n)
{
DWORD end=base+SIFT_BLOCK_LEN-1;
for (i=0;i<SIFT_BLOCK_LEN;i++)
{
g_facArray[i].n=base+i;
g_facArray[i].digSum=0;
}
i=0;
while (primeTab[i]*primeTab[i]<=end && i<primeCount)
{
DWORD k,m;
k= base/primeTab[i];
if (k<primeTab[i])
k=primeTab[i];
m= k*primeTab[i];
if (m<base)
m+= primeTab[i];
while (m <=end)
{
while ( g_facArray[m-base].n % primeTab[i] ==0 )
{
g_facArray[m-base].n /= primeTab[i];
g_facArray[m-base].digSum += getDigSum2(primeTab[i]);
}
m+= primeTab[i];
}
i++;
}
for (i=0;i<SIFT_BLOCK_LEN && i+base<=n;i++)
{
DWORD x=i+base;
WORD s1= getDigSum2(x);
WORD s2= g_facArray[i].digSum;
if ( g_facArray[i].n !=1 && g_facArray[i].n != x)
s2+=getDigSum2(g_facArray[i].n);
if (s1==s2 && s1>0)
{
if (fp!=NULL)
{
numBuff[countInBuff++]=x;
if ( countInBuff*sizeof(DWORD)>=sizeof(numBuff) )
{
fwrite(numBuff,sizeof(DWORD),countInBuff,fp);
countInBuff=0;
}
}
count++;
#ifdef _DEBUG
printf("%u\n",x);
#endif
}
}
base+= SIFT_BLOCK_LEN;
}
if (fp!=NULL && countInBuff!=0 )
{
fwrite(numBuff,sizeof(DWORD),countInBuff,fp);
}
time2=GetTickCount()-time1;
printf("It take %d ms,total %d smith num within %d\n",time2,count,n);
thisExit:
if (primeTab!=NULL)
{
delete[] primeTab; primeTab=NULL;
}
if (fp!=NULL)
{
fclose(fp); fp=NULL;
}
}

复制代码

#include "stdlib.h"
#include "searchSmithNum.h"
void searchSmithNumber(DWORD n,const char *fileName);
int main(int argc, char* argv[])
{
DWORD n=100000000;
searchSmithNumber(n,NULL); //不输出
//searchSmithNumber(n,"smith.dat"); //输出到2进制文件
//计算10^8以内的smith数,不输出，在 PM1.7 需 23.766秒,若输出到话，时间将轻微增加，为24.109秒
//计算10^9以内的smith数,不输出，在 PM1.7 需 252.578秒
return 0;
}

复制代码

无心人 · 发表于 2008-11-27 10:21:36

你这个没达到帖子目的
要求4个以上连续是smith Number的数字的例子

liangbch · 发表于 2008-11-27 10:24:51

该程序在 PIV2.6G 电脑上运行，计算10^8以内的smith数，需 22.128秒。这个CPU
为Northwood 系列，总线 200M，13倍频，512M L2 chche. 需要调整SIFT_BLOCK_LEN 为 (32*1024)，可获得最好的性能。
如果这个程序在主流CPU上运行，速度将更快。大家可给出测试结果。

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

[讨论] K-Smith numbers

回复 25# 无心人的帖子

[讨论] K-Smith numbers

回复 25# 无心人 的帖子

回复 25# 无心人的帖子