分解质因数，生成素数表

经典数论题目详解: POJ_3604-Professor Ben

http://acm.pku.edu.cn/JudgeOnline/problem?id=3604

这个题目足足花费了我一个下午的时间.经历了数十次TLE,RE,CE后,终于将其摆平.

这是一道关于素数的经典好题,之所以这样说是因为它对数学推理和算法复杂度的要求都非常高.先总结一下用到的知识点:

1. 筛法生成素数表(次要)

2. 积性函数和立方公式的推导(核心一)

3. 5万组数据和5百万复杂度的限制下快速分解质因数(核心之二)

 

 题意是给定一个N，求N的所有约数的约数个数的立方之和.(仔细推敲这句话).

首先想到的算法是依次求出从1到5000000每个数约数的个数,然后再利用递推求出结果,但是很不幸严重超时了.

POJ的牛人给出了一个终极公式:

 

其中:

 

 

下面进行的数学推理:

上面公式的意思是将N分解质因数并写成指数的形式,那么最后的结果为这些指数的求和,平方,连乘运算.

1.我们首先假设f(x)代表x的约数的个数,一个数x 的素因数分解为 PI( pi^ki ) . 那么f(x) = PI( 1+ki ).    (PI代表连乘符号)

2.假设题目所求结果为g(x), 则如果x= pi^ki,

则根据题目中的定义易得 g( pi^ki) =f(pi^0)+ f(pi^1) + f(pi^2) +…f(pi^ki) =1^3 +2^3+3^3….(ki+1)^3

这里有一个立方和公式:

1^3 + 2^3 + …… n^3 = [n (n+1) / 2]^2

故上式继续等于 [(1+ki)*(2+ki)]^2

所以我们得到 g( pi^ki) = [(1+ki)*(2+ki)]^2;

现在对于任意一个 x = PI( pi^ki ) ,

易证g(n)为积性函数，即g(nm)=g(n)*g(m)

所以g(x) = g(PI( pi^ki )) = PI( [(1+ki)*(2+ki)]^2);

这就推导出了最后的公式.

 

推出公式后,我们发现要求出结果,必须将5000000以内的数快速分解质因数并且写成指数的形式.

这个也是一个很费神的工作.首先我们生成一个素数表,范围在sqrt(5000000) = 2237以内即可.这个是不超时的关键,因为如果一个数是质数,只需尝试到开方范围内即可,质因数分解时同理.

然后从最小的素数开始去试除一个数,直到全部分解为止,如果超出sqrt(5000000) 仍未分解出因子,说明这个数一定是质数.

 

到此为止这个题也就差不多了,确实是一道很经典的题目,需要认真的推敲.做了一下午,现在头晕眼花了,大伙加油!

 

#include <stdio.h>
#include <math.h>
#include <malloc.h>
static int p[5000];
int prime(int a[],int n) //素数筛法
{
    int i,j,k,x,num,*b;
    n++;
    n/=2;
    b=(int *)malloc(sizeof(int)*(n+1)*2);
    a[0]=2;a[1]=3;num=2;
    for(i=1;i<=2*n;i++)
        b[i]=0;
    for(i=3;i<=n;i+=3)
        for(j=0;j<2;j++)
            {
            x=2*(i+j)-1;
            while(b[x]==0)
                {
                a[num++]=x;
                for(k=x;k<=2*n;k+=x)
                    b[k]=1;
                }
            }
    return num;
}

 

int cal(int n)
{
       int i, ans, back_n, temp;

        back_n = n;

        ans = 1;

        for( i=0; i<=2239 && p[i]*p[i]<=n; ++i ){

                int k;

                if( n%p[i]==0 ){

                        k=0;

                        while( n%p[i]==0 ) k++, n/=p[i];

                        temp =( k+2)*(k+1)/2;

                        ans *= (temp*temp );

                }

        }

        if( n>1 ) ans *=9;
  return ans;
}

int main()
{
 int n;
 int i,j,k;
 int num;

 int p_num = prime(p,3000);

 freopen(”input.txt”,”r”,stdin);
 freopen(”output.txt”,”w”,stdout);

 scanf(”%d”, &n);
 for (i=0; i<n; i++)
 {
  scanf(”%d”, &num);
  printf(”%dn”, cal(num));
 }
 return 0;
}