SPOJ PGCD - Primes in GCD Table (好题! 莫比乌斯反演＋分块求和优化)

PGCD - Primes in GCD Table

Johnny has created a table which encodes the results of some operation -- a function of two arguments. But instead of a boring multiplication table of the sort you learn by heart at prep-school, he has created a GCD (greatest
common divisor) table! So he now has a table (of height a and width
b), indexed from (1,1) to (a,b), and with the value of field (i,j) equal to gcd(i,j). He wants to know how many times he has used prime numbers when writing the table.

Input

First, t ≤ 10, the number of test cases. Each test case consists of two integers, 1 ≤a,b < 10⁷.

Output

For each test case write one number - the number of prime numbers Johnny wrote in that test case.

Example

Input:
2
10 10
100 100

Output:
30
2791

题目链接：http://www.spoj.com/problems/PGCD/

题目大意：1 <= x <= n，1 <= y <= m，求gcd(x, y)为素数的对数

题目分析：和前一题比只是多了一个上界，难度立马变大

学习了这篇博客才解决了这题http://www.cnblogs.com/iwtwiioi/p/4132095.html

为满足且和的的对数

为满足且和的的对数

那么，很显然，反演后得到

因为题目要求是为质数，那么我们枚举每一个质数，然后得到

∑p是质数n∑1<=d<=n/pμ(d)×?n/pd?×?m/pd?

设T=pd，那么问题可以转换为：

∑p是质数n∑1<=d<=n/pμ(d)×?nT?×?mT?

将问题转换为枚举T得：

∑T=1n∑p≤n且p|T且p是质数μ(Tp)×?nT?×?mT?

化简得

∑T=1n?nT?×?mT?∑p≤n且p|T且p是质数μ(Tp)

设

g[x]=∑p≤n且p|x且p是质数μ(xp)

那么原式变成

∑T=1n?nT?×?mT?×g[T]

那么我们只需要考虑如何计算g[x]即可！而且如果p不是质数或者p?T，那么g[x]=0

我们发现g[x]似乎可以在线性筛的时候预处理出？x=k×p,p为质数的情况，可得：

g[kp]={μ(k)μ(k)?g[k]当p|k时当p?k时

首先根据定义，此时

g[kp]=∑p′≤n且p′|kp且p′是质数μ(kpp′)

首先考虑p|k时，有kp质因子p的指数>=2

1、当p′=p，那么约掉后还剩下μ(k)

2、当p′≠p，那么p的质数>=2，根据莫比乌斯函数的定义，为0

因此式1+式2=μ(k)

考虑p?k时，有kp质因子p的指数=1

1、当p′=p，那么约掉后同样是μ(k)

2、当p′≠p，那么因为μ是积性函数，且p与kp′互质，那么得到μ(p)μ(kp′)，然后提出和式。因为μ(p)=?1，而和式内的和恰好就是g[k]的定义，因此这种情况的个数为?g[k]

因此式1+式2=μ(k)?g[k]

最后分块求和然后乘起来

#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;
int const MAX = 1e7 + 5;
int mob[MAX], p[MAX], g[MAX], sum[MAX];
bool prime[MAX];

int Mobius()
{
    mob[1] = 1;
    int pnum = 0;
    for(int i = 2; i < MAX; i++)
    {
        if(!prime[i])
        {
            p[pnum ++] = i;
            mob[i] = -1;
            g[i] = 1;
        }
        for(int j = 0; j < pnum && i * p[j] < MAX; j++)
        {
            prime[i * p[j]] = true;
            if(i % p[j] == 0)
            {
                mob[i * p[j]] = 0;
                g[i * p[j]] = mob[i];
                break;
            }
            mob[i * p[j]] = -mob[i];
            g[i * p[j]] = mob[i] - g[i];
        }
        sum[i] = sum[i - 1] + g[i];
    }
}

ll cal(int l, int r)
{
    ll ans = 0;
    if(l > r)
        swap(l, r);
    for(int i = 1, last = 0; i <= l; i = last + 1)
    {
        last = min(l / (l / i), r / (r / i));
        ans += (ll) (l / i) * (r / i) * (sum[last] - sum[i - 1]);
    }
    return ans;
}

int main()
{
    Mobius();
    int T;
    scanf("%d", &T);
    while(T--)
    {
        int l, r;
        scanf("%d %d", &l, &r);
        printf("%lld\n", cal(l, r));
    }
}

版权声明：本文为博主原创文章，未经博主允许不得转载。