bzoj2820 YY的GCD

Description 神犇YY虐完数论后给傻×kAc出了一题给定N, M,求1<=x<=N, 1<=y<=M且gcd(x,
y)为质数的(x, y)有多少对kAc这种 傻×必然不会了，于是向你来请教……多组输入 Input 第一行一个整数T
表述数据组数接下来T行，每行两个正整数，表示N, M Output T行，每行一个整数表示第i组数据的结果

$$\begin{align} & \sum_p \sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=p] \\ = & \sum_p \sum_{i=1}^{\lfloor\frac np\rfloor}\sum_{j=1}^{\lfloor\frac mp\rfloor}e(\gcd(i,j)) \\ = & \sum_p \sum_{i=1}^{\lfloor\frac np\rfloor}\sum_{j=1}^{\lfloor\frac mp\rfloor} \sum_{x|\gcd(i,j)} \mu(x) \\ = & \sum_p \sum_{i=1}^{\lfloor\frac np\rfloor}\sum_{j=1}^{\lfloor\frac mp\rfloor} \sum_{x|i,x|j} \mu(x) \\ = & \sum_p \sum_{x=1}^{\min(\lfloor\frac np\rfloor,\lfloor\frac mp\rfloor)} \mu(x) \lfloor\frac{\lfloor\frac np\rfloor} x\rfloor \lfloor\frac{\lfloor\frac np\rfloor} x\rfloor \\ = & \sum_p \sum_{x=1}^{\min(\lfloor\frac np\rfloor,\lfloor\frac mp\rfloor)} \mu(x) \lfloor\frac n {px}\rfloor \lfloor\frac m {px}\rfloor \\ = & \sum_{k=1}^{\min(n,m)} \sum_{p是k的质因子} \mu(\frac k p) \lfloor\frac n k\rfloor \lfloor\frac m k\rfloor \end{align}$$
注意到
$$f(k)=\sum_{p是k的质因子} \mu(\frac k p)$$
是一个只跟 $k$有关的式子，如果我们能线性求出$f(k)$，就可以下底函数分块 $O(n+q\sqrt n)$解决问题。
但是要注意 $f(k)$并不是欧拉函数，甚至不是积性函数。但是利用 $\mu(x)$的积性，我们仍然可以线性筛求出 $f(k)$。
具体做法是记下每个数最小质因子的幂，根据线性筛的性质可以知道那个数一定是 $p$，这样记$i=p^mq$，有
$$\begin{align} f(p^mq*p) & = \sum_{x是p^{m+1}q的质因子}\mu(\frac {p^{m+1}q} x) \\ & = \sum_{x是q的质因子}\mu (\frac{p^{m+1}q} x) +\mu(i)\\ & = \sum_{x是q的质因子}\mu (\frac q x)\mu(p^{m+1}) +\mu(i)\\ & =f(q)\mu(p^{m+1})+\mu(i) \\ \end {align}$$
至此，问题解决。

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define LL long long
const int maxn=10000000;
int prm[10000010],mnpow[10000010],mu[10000010],f[10000010],tot;
LL sf[10000010];
void init()
{
    int x;
    for (int i=2;i<=maxn;i++)
    {
        if (!mnpow[i])
        {
            prm[++tot]=i;
            mnpow[i]=i;
            mu[i]=-1;
            f[i]=1;
        }
        for (int j=1;j<=tot&&(LL)i*prm[j]<=maxn;j++)
        {
            x=i*prm[j];
            if (i%prm[j])
            {
                mnpow[x]=prm[j];
                mu[x]=-mu[i];
                f[x]=mu[i]-f[i];
            }
            else
            {
                mnpow[x]=mnpow[i]*prm[j];
                mu[x]=0;
                f[x]=f[i/mnpow[i]]*mu[mnpow[x]]+mu[i]; 
                break;
            }
        }
        sf[i]=sf[i-1]+f[i];
    }
}
LL solve(int n,int m)
{
    LL ret=0;
    int x;
    for (int i=1;i<=n;i=x+1)
    {
        x=min(n/(n/i),m/(m/i));
        ret+=(sf[x]-sf[i-1])*(n/i)*(m/i);
    }
    return ret;
}
int main()
{
    int T,n,m;
    init();
    scanf("%d",&T);
    while (T--)
    {
        scanf("%d%d",&n,&m);
        if (n>m) swap(n,m);
        printf("%lld\n",solve(n,m));
    }
}