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));
    }
}