[BZOJ2820]-YY的GCD-画柿子

说在前面

这题真有意思，突然感觉数论也有套路可循
然而me现在才会做= =…是不是已经没救了

---

题目

BZOJ2820传送门

题面

∑i=1N∑j=1M [gcd(i,j)是质数] $\underset{i=1}{\overset{N}{\sum }}\underset{j=1}{\overset{M}{\sum }}[gcd(i,j)是质数]$
不超过10000组数据

输入输出格式

输入格式：
第一行一个整数T，表示数据组数
接下来T行，每行两个整数N，M表示一组询问

输出格式：
每组数据一行一个整数，表示答案

---

解法

（下面均有 M $M$小于 N $N$）
首先拿到这个柿子，肯定是惯例的化一化
然而发现其中「 gcd(i,j)是质数 $gcd(i,j)是质数$ 」这个条件，好像没有办法表示出来。于是不妨假设我们已经知道了所有小于 M $M$的质数，假设有 k $k$个，质数表示为p，于是有：

∑i=1N∑j=1M[gcd(i,j)为质数]⟹∑i=1N∑j=1M∑d=1k[gcd(i,j)=p(d)]⟹ $$\underset{i=1}{\overset{N}{\sum }}\underset{j=1}{\overset{M}{\sum }}[gcd(i,j)为质数]⟹\underset{i=1}{\overset{N}{\sum }}\underset{j=1}{\overset{M}{\sum }}\underset{d=1}{\overset{k}{\sum }}[gcd(i,j)=p(d)]⟹$$

把d提到前面来，然后套用mobius函数的性质，得到： $$把d提到前面来，然后套用mobius函数的性质，得到：$$

∑d=1k∑i=1Np(d)∑j=1Mp(d)[gcd(i,j)=1]⟹∑d=1k∑i=1Np(d)∑j=1Mp(d)∑t∣i,jμ(t)⟹∑d=1k∑t=1Mp(d)μ(t)∑t∣iNp(d)∑t∣jMp(d)⟹∑d=1k∑t=1Mp(d)μ(t)⌊Np(d)t⌋⌊Mp(d)t⌋⟹ $$\begin{gathered} \underset{d=1}{\overset{k}{\sum }}\underset{i=1}{\overset{\frac{N}{p(d)}}{\sum }}\underset{j=1}{\overset{\frac{M}{p(d)}}{\sum }}[gcd(i,j)=1]⟹\underset{d=1}{\overset{k}{\sum }}\underset{i=1}{\overset{\frac{N}{p(d)}}{\sum }}\underset{j=1}{\overset{\frac{M}{p(d)}}{\sum }}\underset{t\mid i,j}{\sum }\mu (t)⟹ \\ \underset{d=1}{\overset{k}{\sum }}\underset{t=1}{\overset{\frac{M}{p(d)}}{\sum }}\mu (t)\underset{t\mid i}{\overset{\frac{N}{p(d)}}{\sum }}\underset{t\mid j}{\overset{\frac{M}{p(d)}}{\sum }}⟹\underset{d=1}{\overset{k}{\sum }}\underset{t=1}{\overset{\frac{M}{p(d)}}{\sum }}\mu (t)\lfloor \frac{N}{p(d)t}\rfloor \lfloor \frac{M}{p(d)t}\rfloor ⟹ \end{gathered}$$

然后是经典操作：把p(d)t用T来代替，先确定T，再确定p(d) $$然后是经典操作：把p(d)t用T来代替，先确定T，再确定p(d)$$

∑T=1M⌊NT⌋⌊MT⌋∑p(d)∣Tμ(Tp(d)) $$\underset{T=1}{\overset{M}{\sum }}\lfloor \frac{N}{T}\rfloor \lfloor \frac{M}{T}\rfloor \underset{p(d)\mid T}{\sum }\mu \left(\frac{T}{p(d)}\right)$$

到了这一步之后，发现有分块的标志「向下取整符号」，如果我们可以快速的得出 ∑p(d)∣Tμ(Tp(d)) $\underset{p(d)\mid T}{\sum }\mu \left(\frac{T}{p(d)}\right)$的话，那么一个询问就可以在根号时间内被解决了。而且这样做，可以把之前为了方便而引入的 p(d) $p(d)$从最终的式子中移出去

那么设 F(n)=∑p(d)∣nμ(np(d)) $F(n)=\underset{p(d)\mid n}{\sum }\mu \left(\frac{n}{p(d)}\right)$，我们考虑一下如何快速的求出 F(n) $F(n)$
首先我们从含义上去解读这个式子，即「一个整数除掉一个质因数后 μ $\mu$的和」。那么如果 n $n$有两个平方因子 或者 n $n$ 有一个高次方因子，那么 F(n) $F(n)$就为 0 $0$了，据此，我们简单推导一下。
考虑用线性筛的方式，由 F(n) $F(n)$和一个质数 p $p$推出 F(n∗p) $F(n\ast p)$ （即使 F $F$并没有积性，但是我们发现它和 μ $\mu$神似）

• n是p的倍数时
  
  • 当 μ(n)=0 $\mu (n)=0$时，说明其中已经有一个平方因子， p $p$又是一个平方因子，于是 F(n∗p)=0 $F(n\ast p)=0$
  • 不然， p $p$是第一个平方因子，只有 p(d)=p $p(d)=p$的时候， μ(n∗pp(d)) $\mu \left(\frac{n\ast p}{p(d)}\right)$才有值，即 F(n∗p)=μ(n) $F(n\ast p)=\mu (n)$
• n不是p的倍数时，相当于多了一个因子，原来的 μ $\mu$都取反，而 p(d)=p $p(d)=p$时就是 μ(n) $\mu (n)$，即 F(n∗p)=−F(n)+μ(n) $F(n\ast p)=-F(n)+\mu (n)$

那么直接O(n)的求出F()，再求个前缀和，对于每个询问就可以根号出解了

---

下面是自带大常数的代码

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;

bool isnotp[10000005] ;
short mu[10000005] , F[10000005] ;
int T , p[3000000] , pcnt , sF[10000005] , maxM ;
struct Queries{
int N , M ;
}q[10005] ;

void preWork(){
    mu[1] = 1 ;
for( register int i = 2 , j ; i <= maxM ; i ++ ){
if( !isnotp[i] ){
            p[++pcnt] = i ;
            F[i] = 1 , mu[i] = -1 ;
        }
for( j = 1 ; j <= pcnt && i * p[j] <= maxM ; j ++ ){
            isnotp[ i*p[j] ] = true ;
if( i % p[j] == 0 ){
                F[ i*p[j] ] = ( mu[i] == 0 ? 0 : mu[i] ) ;
                mu[ i*p[j] ] = 0 ;
break ;
            }
            F[ i*p[j] ] = -F[i] + mu[i] ;
            mu[ i*p[j] ] = -mu[i] ;
        }
    }
for( register int i = 1 ; i <= maxM ; i ++ )
        sF[i] = sF[i-1] + F[i] ;
}

void solve(){
    register int st , ed ;
for( int k = 1 , N , M ; k <= T ; k ++ ){
        long long ans = 0 ;
        N = q[k].N , M = q[k].M ;
for( st = 1 ; st <= M ; st = ed + 1 ){
            ed = min( N / ( N / st ) , M / ( M / st ) ) ;
            ans += 1LL * ( N / st ) * ( M / st ) * ( sF[ed] - sF[st-1] ) ;
        } printf( "%lld\n" , ans ) ;
    }
}

int main(){
    scanf( "%d" , &T ) ;
for( int i = 1 ; i <= T ; i ++ ){
        scanf( "%d%d" , &q[i].N , &q[i].M ) ; 
if( q[i].N < q[i].M ) swap( q[i].N , q[i].M ) ;
        maxM = max( maxM , q[i].M ) ;
    }
    preWork() ;
    solve() ;
}