liangbch 发表于 2011-11-9 21:20:43

在吃晚饭的时候,回顾了一遍我写的代码,发现一个大bug,待会儿给出修改后的代码

liangbch 发表于 2011-11-9 21:56:29

下面是修正后的代码,速度基本和3楼同// URL: http://bbs.emath.ac.cn/redirect.php?tid=3786&goto=lastpost#lastpost
#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef unsigned long DWORD;
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


static LARGE_INTEGER freq;

static BOOL initFreq()
{
        BOOL ret;
        if ( !QueryPerformanceFrequency( &freq) )
        {        ret=FALSE;        }
        else
        {        ret=TRUE;        }
        return ret;
}

double currTime() //使用高精度计时器
{       
        LARGE_INTEGER performanceCount;
        BOOL result;
        double time=0.0;
        BOOL bRet=TRUE;

        if (freq.QuadPart==0)
        {
                bRet=initFreq();
        }
       
        if (bRet)
        {
                result=QueryPerformanceCounter(&performanceCount );
                time= performanceCount.HighPart * 4294967296.0 + performanceCount.LowPart;
                time=time / (   freq.HighPart * 4294967296.0 + freq.LowPart);
        }
        return time;
}

#define MAX_PRIME 302

struct FAC_ARR_ST
{
        int fc;
        struct
        {
                DWORD p;
                int   e;
        }arr; //2^32以内的数最大9个素因子
};


struct FAC_ARR_ST g_facTab;

void add_to_facTab(DWORD p)
{
        int i;
        int bfind=0;
        for (i=0;i<g_facTab.fc;i++)
        {
                if ( g_facTab.arr[ i ].p==p)
                {
                        g_facTab.arr[ i ].e++;
                        bfind=1;
                        break;
                }
        }

        if ( !bfind )
        {
                g_facTab.arr.p=p;
                g_facTab.arr.e=1;
                g_facTab.fc++;
        }
}

//**********************************************
// n 不能大于L2 cache,否则效能很低
DWORD Sift32(DWORD n,DWORD *primeTab)
// 对0到n 之间(不包括n)的数进行筛选,返回质数的个数
//**********************************************
{
        BYTE *pBuff= new BYTE;
        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; //素数个数

        if ( primeTab==NULL)
        {
                for (i=2;i<n;i++)
                {
                        if (pBuff[ i ])
                                pn++;
                }
        }
        else
        {
                for (i=2;i<n;i++)
                {
                        if (pBuff[ i ])
                                primeTab=i;
                }
        }
        delete[] pBuff;
        return pn;
}


void check_all(DWORD max_n)
{
        DWORD *primeTab=NULL;
        DWORD primeCount;
        DWORD i,j,k;
        int f_c;   //因子数
        int max_fc;//最大因子数

        //创建一个素数表
        primeCount=Sift32(MAX_PRIME,NULL);//得到MAX_PRIME以内的素数个数
        primeTab=(DWORD *)malloc(sizeof(DWORD)* (primeCount+1));
        Sift32(MAX_PRIME,primeTab);
        primeTab=2147483647; //设置一个栅栏

       
        for (max_fc=0,i=2;i<=max_n;i++)
        {
                DWORD na;
                if (i & 1) //i是奇数
                {    na=i; na=(i+1)/2;}
                else
                {   na=i/2; na=i+1;}
               
                for (g_facTab.fc=0,j=0;j<2;j++)
                {
                        DWORD n,r;
                       
                        r=n=na;
                        // primeTab=2147483647,保证循环能够终止
                       
                        for (k=0; primeTab * primeTab <= n ; )
                        {
                                DWORD prime=primeTab;
                                if ( r % prime==0)
                                {
                                        while ( r % prime==0)
                                        {
                                                add_to_facTab(prime);//增加素因子到素因子表 g_facTab
                                                r/=prime;
                                        }
                                }
                                k++;
                                if ( primeTab > r)
                                        break;
                        }

                        if (r >1)
                                add_to_facTab(r);//增加素因子到素因子表 g_facTab
                }

                for (f_c=1,j=0;j<g_facTab.fc;j++)
                        f_c *= (g_facTab.arr.e +1);

                //printf("%u*(%u+1)/2 with %u factor\n",i,i,f_c);
                if (f_c >max_fc)
                {
                        printf("%u*(%u+1)/2 with %u factor\n",i,i,f_c);
                        max_fc=f_c;
                }
        }
        free(primeTab);
}

int main(int argc, char* argv[])
{
        double t=currTime();
        //check_all(92681);
        check_all(65535);
        t=currTime()-t;
        printf("it take %.6f ms\n",t*1000);
        return 0;
}

liexi20101117 发表于 2011-11-12 22:14:59

版主的代码需要的时间啃下哟,学习中

shingyuan 发表于 2011-11-29 23:07:36

非常感谢各位朋友的帮助,这些代码需要好好消化学习
页: 1 [2]
查看完整版本: 求因子个数的问题