首页 > 技术文章 > SPOJ PGCD - Primes in GCD Table (好题! 莫比乌斯反演+分块求和优化)

jzdwajue 2017-08-13 14:09 原文


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 < 107.

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


为满足的对数

为满足的对数

 那么,非常显然,反演后得到

由于题目要求是为质数,那么我们枚举每个质数,然后得到

pn1<=d<=n/pμ(d)×n/pd×m/pd

T=pd。那么问题能够转换为:

pn1<=d<=n/pμ(d)×nT×mT

将问题转换为枚举T得:

T=1npnp|Tpμ(Tp)×nT×mT

化简得

T=1nnT×mTpnp|Tpμ(Tp)

g[x]=pnp|xpμ(xp)

那么原式变成

T=1nnT×mT×g[T]

那么我们仅仅须要考虑怎样计算g[x]就可以。并且假设ppT。那么g[x]=0

我们发现g[x]似乎能够在线性筛的时候预处理出?x=k×p,p的情况,可得:

g[kp]={μ(k)μ(k)g[k]p|kpk

首先依据定义。此时

g[kp]=pnp|kppμ(kpp)

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

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

2、当pp。那么p的质数>=2。依据莫比乌斯函数的定义。为0

因此1+2=μ(k)

考虑pk时,有kp质因子p的指数=1

1、当p=p。那么约掉后相同是μ(k)

2、当pp,那么由于μ是积性函数,且pkp互质,那么得到μ(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));
    }
}

推荐阅读