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

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

:lol

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

回复 25# 无心人 的帖子

使用晒法,我在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;
        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)
                {
                        primeTab=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;
FAC_ST        g_facArray;

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);

        DWORD s=0;
        if (n<10000)
                return g_digSumTab;
        else if (n<100000000)
        {
                DWORD t=n % 10000;
                n/=10000;
                return g_digSumTab+g_digSumTab;
        }
        else
        {
                DWORD t0,t1;
                t0=n % 10000;
                n/=10000;
                t1=n % 10000;
                n/=10000;
                return g_digSumTab+g_digSumTab+g_digSumTab;
        }
}

void initTable(BYTE *tab,DWORD n)
{
        DWORD i;
        assert(n==10 || n==100 || n==1000 || n==10000);
        for (i=0;i<n;i++)
        {
                tab=(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;
        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;
        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.n=base+i;
                        g_facArray.digSum=0;
                }
               
                i=0;
                while (primeTab*primeTab<=end && i<primeCount)
                {
                        DWORD k,m;
                       
                        k= base/primeTab;
                        if (k<primeTab)
                                k=primeTab;
                        m= k*primeTab;
                        if (m<base)
                                m+= primeTab;
                       
                        while (m <=end)
                        {
                                while ( g_facArray.n% primeTab ==0 )
                                {
                                        g_facArray.n /= primeTab;
                                        g_facArray.digSum += getDigSum2(primeTab);
                                }
                                m+= primeTab;
                        }
                        i++;
                }
               
                for (i=0;i<SIFT_BLOCK_LEN && i+base<=n;i++)
                {
                        DWORD x=i+base;
                        WORD s1= getDigSum2(x);
                        WORD s2= g_facArray.digSum;
                       
                        if ( g_facArray.n !=1 && g_facArray.n != x)
                                s2+=getDigSum2(g_facArray.n);
                        if (s1==s2 && s1>0)
                        {
                                if (fp!=NULL)
                                {
                                        numBuff=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

:lol

你这个没达到帖子目的
要求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上运行,速度将更快。大家可给出测试结果。
页: 1 2 [3] 4 5 6 7
查看完整版本: K-Smith numbers