YY的GCD【luoguP2257】


题目大意

有至多(10000)组询问,问(1 < i leqslant N leqslant 10000000, 1 < j leqslant M leqslant 10000000),并且(gcd(i, j))为质数的有多少对。

解题思路

为了方便描述,我们定义([]),当([])中表达式为真时为(1),否则为(0)。同时定义(Prime)为素数集合。
下面的讨论中,我们不妨设(N leqslant M)
我们设

[f(d)=sum_{i=1}^Msum_{j=1}^N[gcd( i, j ) in Prime]\ F(n)=sum_{n|d}^Nf(d)=lfloorfrac{M}{n} floorlfloorfrac{N}{n} floor ]

即,(f(d))是当(gcd=d)时的答案数,(F(n))是当(gcd)(n)的倍数时的答案数。
我们发现,求(F(n))十分的方便,于是我们考虑能否通过(F(n))(f(n))表述出来。

由莫比乌斯反演,得

[f(n)=sum_{n|d}mu(frac{d}{n})F(d) ]

其中(mu)是莫比乌斯函数。

那么答案就可以表示为

[egin{aligned} Ans & = sum_{nin Prime}^Nf(n)\ & = sum_{nin Prime}^Nsum_{n|d}mu(frac{d}{n})F(d)\ & = sum_{nin Prime}^Nsum_{n|d}mu(frac{d}{n})lfloorfrac{N}{d} floorlfloorfrac{M}{d} floor\ & = sum_{d}^Nsum_{n|d,nin Prime}mu(frac{d}{n})lfloorfrac{N}{d} floorlfloorfrac{M}{d} floor\ & = sum_{d}^Nlfloorfrac{N}{d} floorlfloorfrac{M}{d} floorsum_{n|d, nin Prime} mu(frac{d}{n}) end{aligned} ]

通过稍微修改线筛,我们可以与处理出(mu),然后可以预处理出所有的(sum_{n|d,nin Prime}mu(frac{d}{n}))。最后再整除分块统计答案就可以了。

参考程序

程序中,mu即为(mu)(Sum)为前缀和。

#include <bits/stdc++.h>
using namespace std;

const int MaxN = 10000010;
int Mu[ MaxN ], Vis[ MaxN ];
long long Sum[ MaxN ];
int Num, Prime[ 1000010 ];

void Init() {
    Mu[ 1 ] = 1;
    for( int i = 2; i <= MaxN; ++i ) {
        if( !Vis[ i ] ) Prime[ ++Num ] = i, Mu[ i ] = -1;
        for( int j = 1; j <= Num && ( long long ) i * Prime[ j ] <= ( long long ) MaxN; ++j ) {
            Vis[ i * Prime[ j ] ] = 1;
            if( i % Prime[ j ] == 0 ) break;
            Mu[ i * Prime[ j ] ] = - Mu[ i ];
        }
    }
    for( int i = 1; i <= MaxN; ++i ) 
        for( int j = 1; j <= Num && ( long long ) i * Prime[ j ] <= ( long long ) MaxN; ++j )
            Sum[ i * Prime[ j ] ] += Mu[ i ];
    for( int i = 2; i <= MaxN; ++i ) Sum[ i ] += Sum[ i - 1 ];
    return;
}

void Work() {
    int N, M; 
    scanf( "%d%d", &N, &M );
    if( N > M ) swap( N, M );
    long long Ans = 0;
    for( int x = 1, y; x <= N; x = y + 1 ) {
        y = min( N / ( N / x ), M / ( M / x ) );
        Ans += 1LL * ( N / x ) * ( M / x ) * ( Sum[ y ] - Sum[ x - 1 ] );
    }
    printf( "%lld
", Ans );
    return;
}

int main() {
    Init();
    int T; scanf( "%d", &T );
    for( ; T; --T ) Work();
    return 0;
}
原文地址:https://www.cnblogs.com/chy-2003/p/10138031.html