首页 > ACM题库 > HDU-杭电 > HDU 4746-Mophues-数论-[解题报告]HOJ
2015
09-17

HDU 4746-Mophues-数论-[解题报告]HOJ

Mophues

问题描述 :

As we know, any positive integer C ( C >= 2 ) can be written as the multiply of some prime numbers:
    C = p1×p2× p3× … × pk
which p1, p2 … pk are all prime numbers.For example, if C = 24, then:
    24 = 2 × 2 × 2 × 3
    here, p1 = p2 = p3 = 2, p4 = 3, k = 4

Given two integers P and C. if k<=P( k is the number of C’s prime factors), we call C a lucky number of P.

Now, XXX needs to count the number of pairs (a, b), which 1<=a<=n , 1<=b<=m, and gcd(a,b) is a lucky number of a given P ( "gcd" means "greatest common divisor").

Please note that we define 1 as lucky number of any non-negative integers because 1 has no prime factor.

输入:

The first line of input is an integer Q meaning that there are Q test cases.
Then Q lines follow, each line is a test case and each test case contains three non-negative numbers: n, m and P (n, m, P <= 5×105. Q <=5000).

输出:

The first line of input is an integer Q meaning that there are Q test cases.
Then Q lines follow, each line is a test case and each test case contains three non-negative numbers: n, m and P (n, m, P <= 5×105. Q <=5000).

样例输入:

2
10 10 0
10 10 1

样例输出:

63
93

题目:http://acm.hdu.edu.cn/showproblem.php?pid=4746

 

题意:给出n, m, p,求有多少对a, b满足gcd(a, b)的素因子个数<=p,(其中1<=a<=n, 1<=b<=m)

分析:设A(d):gcd(a, b)=d的有多少种

     设B(j): gcd(a, b)是j的倍数的有多少种,易知B(j) = (n/j)*(m/j)

     则由容斥原理得:(注:不同行的μ是不相同的,μ为莫比乌斯函数)

     A(1) = μ(1)*B(1) + μ(2)*B(2) + μ(3)*B(3) + … + μ(p1*p2…)*B(p1*p2…)

     A(2) = μ(1)*B(1*2) + μ(2)*B(2*2) + μ(3)*B(3*2) + … + μ(p1*p2..)*B(p1*p2..*2)

     …

     A(d) = μ(1)*B(1*d) + μ(2)*B(2*d) + μ(3)*B(3*d) + … + μ(p1*p2..)*B(p1*p2..*d)

 

     ans = A(1)+A(2)+…+A(d) = F(1)*B(1) + F(2)*B(2) + … + F(p1*p2..)*B(p1*p2..)

     于是可以枚举公约数i{表示A(i)},利用筛法找出i的倍数j,i对B(j)的贡献系数为:F(j)+=μ(j/i)

     总之,求出B(j)的总贡献系数F(j)即可得答案:F(1)*B(1)+F(2)*B(2)+…+F(n)*B(n)

     上面没有限制gcd的素因子个数,要限制其实不难,给系数加多一维即可:

     F(d)(p)表示:素因子个数<=p时,对B(d)的贡献系数

   

     分块加速思想

     你可以再纸上模拟一下:设d在[i, n/(n/i)]的区间上,则该区间内所有的n/d都是一样的。

 

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
#define LL long long
#define M 500005
#define N 19

//返回n中有多少个x因子
int cal(int n, int x)
{
    int res = 0;
    do
    {
        ++res;
        n /= x;
    }
    while (n % x == 0);
    return res;
}

//备注:分块加速求解需要求前缀和
//F[i][j]: 表示素因子个数<=j条件下的莫比乌斯前缀和:μ(1)+μ(2)+...+μ(i)
int F[M][N];
int num[M];		//num[i]: i中含有多少个素因子
int h[M];		//h[i]: -1表示存在平方因子,否则表示有多少种素因子

//莫比乌斯函数的定义
int mob(int n)
{
    if (h[n] == -1) return 0;	//存在平方因子时,μ(n)=0
    if (h[n] & 1) return -1;	//奇数个不同素数之积,μ(n)=-1
    return 1;					//偶数个不同素数之积,μ(n)=1
}

int main()
{
    int t, n, m, p, i, j;
    //筛法算出num[]以及h[]
    for (i = 2; i < M; i++)
    {
        if (num[i]) continue;
        for (j = i; j < M; j+=i)
        {
            int tp = cal(j, i);
            num[j] += tp;
            if (tp > 1)  	//j中含有多个i,必然存在平方因子
            {
                h[j] = -1;
            }
            else if (h[j] >= 0)
            {
                ++h[j];
            }
        }
    }
    //枚举i作为公因子,对B(j)的贡献值为:mob(j/i)
    for (i = 1; i < M; i++)
    {
        for (j = i; j < M; j+=i)
        {
            F[j][num[i]] += mob(j/i);
        }
    }
    //为了表示素因子数<=j的意义,求j的前缀和
    for (i = 1; i < M; i++)
    {
        for (j = 1; j < N; j++)
        {
            F[i][j] += F[i][j-1];
        }
    }
    //为了分组加速求解,求i的前缀和
    for (i = 1; i < M; i++)
    {
        for (j = 0; j < N; j++)
        {
            F[i][j] += F[i-1][j];
        }
    }
    scanf("%d", &t);
    while (t--)
    {
        scanf("%d%d%d", &n, &m, &p);
        LL ans = 0;
        if (p >= N)
        {
            ans = (LL)n*m;
        }
        else
        {
            if (n > m)
            {
                n ^= m;
                m ^= n;
                n ^= m;
            }
            for (i = 1; i <= n; i = j + 1)
            {
                j = min(n/(n/i), m/(m/i));
                ans += ((LL)F[j][p]-F[i-1][p])*(n/i)*(m/i);
            }
        }
        printf("%I64d\n", ans);
    }
    return 0;
}

 

参考:http://blog.csdn.net/acdreamers/article/details/12871643