[SDOI2015]约数个数和 莫比乌斯反演

～～～题面～～～

题解：

为什么SDOI这么喜欢莫比乌斯反演，，，

首先有一个结论$$d(ij) = sum_{x|i}sum_{y|j}[gcd(x, y) == 1]$$
为什么呢？
首先，可以看做从两个数中分别取一些不重叠的质数的$k_{i}$次方，组成新数的方案数。
那如果有需要重叠的部分怎么办？
可以看做全都在第一个or第二个数中取。
但是一个数的次数不够怎么办呢？
相当于以1为媒介，可以简介统计到这些情况
比如$2^{3} cdot 2^{2}  = 2^{5}$可以看成$1, 2^{3}$ + $1, 2^{2}$,这样就可以做到分别枚举了指数从1到5的情况了；
$$ans = sum_{i = 1}^{n}sum_{j = 1}^{m}sum_{x | i}sum_{y | j}[gcd(x, y) == 1]$$
$$ans = sum_{x = 1}^{n}sum_{y = 1}^{m} lfloor{frac{n}{x}} floor lfloor{frac{m}{y}} floor [gcd(x, y) == 1]$$<---改成枚举x和y，然后统计有多少对i,j中分别含有因数x，y
设$$f(d) = sum_{x = 1}^{n}sum_{y = 1}^{m} lfloor{frac{n}{x}} floor lfloor{frac{m}{y}} floor [gcd(x, y) == d]$$, $$g(x) = sum_{x | d}^{min(n, m)}{f(d)}$$
那么$ans = f(1)$
反演一下得到:$$f(n)  = sum_{n | d}mu(frac{d}{n}) cdot g(d)$$
接下来就是怎么求$g(d)$了。
$$g(d) = sum_{x = 1}^{n}sum_{y = 1}^{m} lfloor{frac{n}{x}} floor lfloor{frac{m}{y}} floor [d | gcd(x, y)]$$
$d | gcd(x, y)$ ---> $d | x + d | y$ ,因为$x | n + y | m$ ---> $d | n + d | m$
所以枚举$d_{i}$就符合条件了，$lfloor{frac{n}{dx}} floor$ ---> $dx$倍数的个数，即有多少i，j的搭配是可以满足$[d|gcd(x, y)]$的
$$g(d) = sum_{x = 1}^{ lfloor{frac{n}{d}} floor }sum_{y = 1}^{ lfloor{frac{m}{d}} floor } lfloor{frac{n}{dx}} floor lfloor{frac{m}{dy}} floor$$
于是现在就是要求$f(1)$
$$f(1) = sum_{d = 1}^{min(n, m)}{mu(frac{d}{1})g(d)}$$
$$= sum_{d = 1}^{min(n, m)}{mu(d)g(d)}$$
$$= sum_{d = 1}^{min(n, m)}{mu(d) cdot sum_{x = 1}^{ lfloor{frac{n}{d}} floor }sum_{y = 1}^{ lfloor{frac{m}{d}} floor } lfloor{frac{n}{dx}} floor lfloor{frac{m}{dy}} floor}$$
这个式子可以整数分块。
设$N = frac{n}{d}$， $M = frac{m}{d}$,那么
$$= sum_{d = 1}^{min(n, m)}{mu(d) cdot sum_{x = 1}^{ N }sum_{y = 1}^{ M } lfloor{frac{N}{x}} floor lfloor{frac{M}{y}} floor}$$
$$= sum_{d = 1}^{min(n, m)}{mu(d) cdot sum_{x = 1}^{ N}{lfloor{frac{N}{x}} floor} sum_{y = 1}^{M}  {lfloor{frac{M}{y}} floor}}$$
而$N = lfloor{frac{n}{d}} floor$， $M = lfloor{frac{m}{d}} floor$有很多相同的区间段,所以可以整数分块求。

#include<bits/stdc++.h>
using namespace std;
#define R register int
#define AC 50100
#define ac 50000
#define LL long long
int n, m, tot, t;
LL ans;
int prime[AC];
LL mu[AC], s[AC];
bool z[AC];

inline int read()
{
    int x = 0;char c = getchar();
    while(c > '9' || c < '0') c = getchar();//error!!!又一次打错了读入优化。。。是||啊
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x;
}

void pre()
{
    int x;
    mu[1] = 1;
    for(R i = 2; i <= ac; i++)
    {
        if(!z[i]) prime[++tot] = i, mu[i] = -1;
        for(R j = 1; j <= tot; j++)
        {
            x = prime[j];
            if(x * i > ac) break;
            z[x * i] = true;
            if(!(i % x)) break;
            mu[i * x] = -mu[i]; 
        }
    }
    for(R i = 1; i <= ac; i++) mu[i] += mu[i-1];    
    int pos;
    for(R i = 1; i <= ac; i++) 
    {
        for(R j = 1; j <= i; j = pos + 1)
        {
            pos = i / (i / j);
            s[i] += (LL)(i / j) * (LL)(pos - j + 1);
        }
    }
//    for(R i = 1; i <= 30; i++) printf("%lld
", s[i]);
}

void work()
{
    t = read();
    while(t--)
    {
        n = read(), m = read(), ans = 0;
        //scanf("%d%d", &n, &m);ans = 0;
        int pos, b = min(n, m);
        for(R i = 1; i <= b; i = pos + 1)
        {
            pos = min(n / (n / i), m / (m / i));
            ans += (mu[pos] - mu[i - 1]) * s[n / i] * s[m / i];
        }
        printf("%lld
", ans);
    }
}

int main()
{
//    freopen("in.in", "r", stdin);
    pre();
    work();
//    fclose(stdin);
    return 0;
}