[BZOJ2820]-YY的GCD-画柿子

说在前面

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

---

题目

BZOJ2820传送门

题面

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

输入输出格式

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

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

---

解法

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

$$\sum\limits_{i=1}^N\sum\limits_{j=1}^M[ gcd(i,j)为质数] \implies \sum\limits_{i=1}^N\sum\limits_{j=1}^M\sum\limits_{d=1}^k[gcd(i,j)=p(d)]\implies$$

$$把d提到前面来，然后套用mobius函数的性质，得到：$$ $$\sum\limits_{d=1}^k\sum\limits_{i=1}^\frac{N}{p(d)}\sum\limits_{j=1}^\frac{M}{p(d)}[gcd(i,j)=1]\implies \sum\limits_{d=1}^k\sum\limits_{i=1}^\frac{N}{p(d)}\sum\limits_{j=1}^\frac{M}{p(d)}\sum\limits_{t\mid i,j}\mu(t)\implies \\ \sum\limits_{d=1}^k\sum\limits_{t=1}^{\frac{M}{p(d)} }\mu(t)\sum\limits_{t\mid i}^{\frac{N}{p(d)}}\sum\limits_{t\mid j}^{\frac{M}{p(d)}} \implies \sum\limits_{d=1}^k\sum\limits_{t=1}^{\frac{M}{p(d)} }\mu(t)\left\lfloor \frac{N}{p(d)t}\right\rfloor \left\lfloor \frac{M}{p(d)t}\right\rfloor \implies$$ $$然后是经典操作：把p(d)t用T来代替，先确定T，再确定p(d)$$ $$\sum\limits_{T=1}^M\left\lfloor\frac{N}{T}\right\rfloor \left\lfloor\frac{M}{T}\right\rfloor \sum\limits_{p(d)\mid T}\mu\left(\frac{T}{p(d)}\right)$$

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

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

• n是p的倍数时
  
  • 当$\mu(n)=0$时，说明其中已经有一个平方因子，$p$又是一个平方因子，于是$F(n*p)=0$
  • 不然，$p$是第一个平方因子，只有$p(d)=p$的时候，$\mu\left(\frac{n*p}{p(d)}\right)$才有值，即$F(n*p)=\mu(n)$
• n不是p的倍数时，相当于多了一个因子，原来的$\mu$都取反，而$p(d)=p$时就是$\mu(n)$，即$F(n*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() ;
}