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

---

Added by:	Yash
Date:	2009-06-12
Time limit:	0.687s
Source limit:	11111B
Memory limit:	1536MB
Cluster:	Cube (Intel Pentium G860 3GHz)
Languages:	All except: ERL JS NODEJS PERL 6 VB.net
Resource:	Codechef

题目链接：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 noprime[MAX];

int Mobius()
{
    mob[1] = 1;
    int pnum = 0;
    for(int i = 2; i < MAX; i++)
    {
        if(!noprime[i])
        {
            p[pnum ++] = i;
            mob[i] = -1;
            g[i] = 1;
        }
        for(int j = 0; j < pnum && i * p[j] < MAX; j++)
        {
            noprime[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));
    }
}