SPOJ PGCD - Primes in GCD Table (好題! 莫比烏斯反演＋分塊求和優化)

_TCgogogo_發表於2015-08-20
GC優化


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


為滿足且和的的對數

為滿足且和的的對數

 那麼，很顯然，反演後得到

因為題目要求是為質數，那麼我們列舉每一個質數，然後得到

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

相關文章