Sdoi2015约数个数和题解莫比乌斯反演

题目描述

• T组数据，求$\Sigma_{i=1}^N\Sigma_{j=1}^Md(ij)$，$d(x)$代表x的约数个数。

• $1\le N, M,T\le 10^5$

题解

• 首先，膜拜一下PoPoQQQ大神及其题解

• 然后，有一个神奇的结论：$\Sigma_{i=1}^N\Sigma_{j=1}^Md(ij)$=$\Sigma_{i=1}^N\Sigma_{j=1}^M[\frac{N}{i}][\frac{M}{j}][\gcd(i,j)=1]$
  要证上式，只需证$d(nm)=\Sigma_{i|n}\Sigma_{j|m}[\gcd(i,j)=1]$，因为上式即为该式的前缀和形式。
  分开考虑每个质数$p$对答案的贡献。设$n=n'p^{k_1},m=m'p^{k_2}$，那$p$对$d(nm)$的贡献就是$k_1+k_2+1$，对等式右边的贡献是有序数对$(i,j)$中的$(p^{k_1},1)(p^{k_1-1},1)(p^{k_1-2},1)\dots(1,1)(1,p)\dots(1,p^{k_2-1})(1,p^{k_2})$也是$k_1+k_2+1$，所以原式得证。

• 再令$f(x)=\Sigma_{i=1}^x\left[\frac{x}{i}\right]$，显然有$f(x)=d(x)$，因为$f(x)$相当于把$i=1\to x$对$d(x)$的贡献求和，就是$d(x)$本身

• 然后可以莫比乌斯反演了：
  

  $$ans=\Sigma_{i=1}^N\Sigma_{j=1}^M\left[\frac{N}{i}\right]\left[\frac{M}{j}\right][\gcd\left(i,j\right)=1]\\=\Sigma_{i=1}^N\Sigma_{j=1}^M\Sigma_{k|\gcd\left(i,j\right)}\mu\left(k\right)\left[\frac{N}{i}\right]\left[\frac{M}{j}\right]\\=\Sigma_{k=1}^N\Sigma_{i=1}^{\left[\frac{N}{k}\right]}\Sigma_{j=1}^{\left[\frac{M}{k}\right]}\mu\left(k\right)\left[\frac{N}{ki}\right]\left[\frac{M}{kj}\right]\\=\Sigma_{k=1}^N\mu\left(k\right)\Sigma_{i=1}^{\left[\frac{N}{k}\right]}\left[\frac{N}{ki}\right]\Sigma_{j=1}^{\left[\frac{M}{k}\right]}\left[\frac{N}{kj}\right]\\=\Sigma_{k=1}^N\mu\left(k\right)f\left(\left[\frac{M}{k}\right]\right)f\left(\left[\frac{M}{k}\right]\right)\\=\Sigma_{k=1}^N\mu(k)d\left(\left[\frac{N}{k}\right]\right)d\left(\left[\frac{M}{k}\right]\right)$$

• 然后用筛法筛出$\mu(i)$和$d(i)$再做就好了。

• 看来还需要锻炼自己的数学功底啊

Code

#include <cstdio>
#include <algorithm>
#define ll long long
#define M 50000
using namespace std;
int mu[M + 5], n, m, pri[M], top, d[M + 5], f[M + 5];
bool mrk[M + 5];
void getmu()
{
    mu[1] = d[1] = 1;
    for(int i = 2; i <= M; ++i)
    {
        if(!mrk[i])
        {
            pri[top++] = i;
            mu[i] = -1;
            d[i] = 2;
            f[i] = 1;
        }
        for(int j = 0; j < top && i * pri[j] <= M; ++j)
        {
            mrk[i * pri[j]] = true;
            if(i % pri[j] == 0)
            {
                mu[i * pri[j]] = 0;
                f[i * pri[j]] = f[i] + 1;
                d[i * pri[j]] = d[i] / (f[i] + 1) * (f[i] + 2);
                break;
            }
            mu[i * pri[j]] = -mu[i];
            d[i * pri[j]] = (d[i] << 1);
            f[i * pri[j]] = 1;
        }
    }
    for(int i = 2; i <= M; ++ i)
    {
        d[i] += d[i - 1];
        mu[i] += mu[i - 1];
    }
}
inline ll work(int n, int m)
{
    if(n > m) swap(n, m);
    ll ans = 0;
    for(int i = 1, j; i <= n; i = j + 1)
    {
        j = min(n / (n / i), m / (m / i));
        ans += (ll)(mu[j] - mu[i - 1]) * (ll)d[n / i] * (ll)d[m / i];
    }
    return ans;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("d.txt", "r", stdin);
    freopen("a.txt", "w", stdout);
#endif
    int T;
    getmu();
    scanf("%d", &T);
    while(T--)
    {
        scanf("%d%d", &n, &m);
        printf("%lld\n", work(n, m));
    }
    return 0;
}