Enzymii的hu测 T3.violetevergarden（反演【约数个数）

版权属于Enzymii，想要引用此题（包括题面）的朋友请联系博主

分析：

题目来源

首先要好好读题，题目有一点问题
不过不要在意这么多细节，实际上题目要求：

$$\sum_{i=1}^{n}d(i^2)$$

d的暴力计算公式
$n=\prod p_i^{k_i}$
$d(n)=\sum (k_i+1)$

ZYXZYXZYX和我说有一个60%的做法，简单反演一下就可以了：

首先，关于$d$函数（约束个数）有一个玄学公式

$$d(nm)=\sum_{i|n}\sum_{j|m}[gcd(i,j)=1]$$

$$\sum_{k=1}^{n}d(k^2)=\sum_{k=1}^{n}\sum_{i|k}\sum_{j|k}[gcd(i,j)=1]$$

$$=\sum_{k=1}^{n}\sum_{i|k}\sum_{j|k}\sum_{d|gcd(i,j)}μ(d)$$

$$=\sum_{k=1}^{n}\sum_{i|k}\sum_{j|k}\sum_{d=1}^{k}\sum_{d|i}\sum_{d|j}μ(d)$$

$$=\sum_{k=1}^{n}\sum_{d=1}^{k}μ(d)\sum_{i|k}\sum_{j|k}\sum_{d|i}\sum_{d|j}1$$

$$=\sum_{k=1}^{n}\sum_{d=1}^{k}μ(d)\sum_{i|k}\sum_{d|i}\sum_{j|k}\sum_{d|j}1$$

显然$\sum_{i|k}\sum_{d|i}\sum_{j|k}\sum_{d|j}1$就是$μ(d)$的计算次数
而通过观察我们可以发现$\sum_{i|k}\sum_{d|i}$和$\sum_{j|k}\sum_{d|j}$实际上是等价的

下面我们考虑如何计算$\sum_{i|k}\sum_{d|i}1$
上式的意义：i是k的约数，d是i的约数
设$i=pd$，则有

$$\sum_{pd|k}1=\sum_{p|{k \over d}}1$$

仔细观察上式，意义即为：${k \over d}$的约数个数和，即$d({k \over d})$

因此原式可化为：

$$\sum_{k=1}^{n}\sum_{d=1}^{k}μ(d)d^2({k \over d})$$

满分做法

我们要用到暴力计算d的公式
设$n=\prod_{i=1}^{k}p_i^{a_i}$，则$n^2=\prod_{i=1}^{k}p_i^{2a_i}$

$$d(n)=\sum_{i=1}^{k}(a_i+1)$$

$$d(n^2)=\sum_{i=1}^{k}(2a_i+1)$$

因为$n=\prod_{i=1}^{k}p_i^{a_i}$，那么$n$的约数中$p_i$的幂范围就是$0-a_i$
因为$n^2=\prod_{i=1}^{k}p_i^{2a_i}$，那么$n^2$的约数中$p_i$的幂范围就是$0-2a_i$
$n$的约数一定是$n^2$的约数，那么我们考虑$n^2$的约数$X$，但是$X$不是$n$的约数
那么$X$一定可以表示为

$$X=n的某个约数*\prod p_i^{a_i}$$
意义即为：X比n的某个约数多乘上若干个 $p_i^{a_i}$
那么到底是多乘上了那些 $p_i^{a_i}$呢？我们可以枚举

设$n$的不同约数有$\omega(n)$个
那么这$\omega(n)$个元素能组成的集合个数：

$$2^{\omega(n)}$$

对于上面的$2^{\omega(n)}$集合，我们都可以让ta与$n$的任何一个约数重组出$n^2$的合法约数

举个例子：
$n=p_1^{a_1}p_2^{a_2}$
$n^2=p_1^{2a_1}p_2^{2a_2}$
$\omega(n)=2$
那么我们能得到的集合为$\{\{0,0\},\{0,1\},\{1,0\},\{1,1\}\}$
以n的一个约数$p_1^{1}p_2^{1}$为基准，上面的集合与约数$p_1^{1}p_2^{1}$能组成的$n^2$的约数分别为：
$p_1^{1}p_2^{1},p_1^{1}p_2^{a_2+1},p_1^{a_1+1}p_2^{1}，p_1^{a_1+1}p_2^{a_2+1}$

考虑如何快速的计算$2^{\omega(n)}$
观察这些集合的特点，集合元素非0即1

举个例子：
$n=p_1^{a_1}p_2^{a_2}p_3^{a_3}$
$\omega(n)=2$
一个集合为$\{0,1,1\}$
显然$\{0,1,1\}$分别对应着$p_1,p_2,p_3$是否选择了
如果我们干脆把$\{0,1,1\}$视为$p_1,p_2,p_3$的指数
得到一个数$p_1^0+p_2^0+p_3$

我们用上述方法得到的的数每个质因子的指数非0即1，也就是说没有平方因子
一提到没有平方因子，就可以想到$μ$函数，实际上：

$$2^{\omega(n)}=\sum_{d|n}μ^2(d)$$

意义：$n$的没有平方因子（$μ!=0$）的约数的计数

回到我们一开始的问题：

$$d(n^2)=\sum_{i=1}^{k}(2a_i+1)$$

之前我们说过：对于$2^{\omega(n)}$个集合，我们都可以让ta与$n$的任何一个约数重组出$n^2$的合法约数

$$d(n^2)=\sum_{i=1}^{k}(2a_i+1)=\sum_{i|n}d(i)*2^{\omega(n)}=\sum_{i|n}d(i)\sum_{j|n}μ^2(j)$$

枚举$μ$：

$$\sum_{i=1}^{n}μ^2(i)\sum_{j=1}^{n \over i}d(j)$$

于是利用整除分块，我们只要知道$μ^2(i)$和$d(i)$的前缀和，预处理前$n^{2 \over 3}$项的情况下，$O(n^{2 \over 3})$求出

前者就是$1-n$的无平方因子数个数，考虑容斥$O(\sqrt n)$

$$\sum_{i=1}^{n}μ^2(i)=\sum_{i=1}^{\sqrt n}μ(i)*{n \over i^2}$$

后者就是套路

$$\sum_{i=1}^{n}d(i)=\sum_{i=1}^{n}\sum_{j=1}^{1}{n \over i}{1 \over j}[gcd(i,j)=1]=\sum_{i=1}^{n}{n \over i}$$

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#define ll long long

using namespace std;

const int N=10000001;
ll n;
int sshu[N],tot=0;
ll mu[N],d[N],sm[N],sd[N],c[N];
bool no[N];

void prepare() {
    mu[1]=1; d[1]=1;
    for (int i=2;i<N;i++) {
        if (!no[i]) {
            sshu[++tot]=i;
            mu[i]=-1;
            d[i]=2; c[i]=1;
        }
        for (int j=1;j<=tot&&sshu[j]*i<N;j++) {
            int v=sshu[j]*i;
            no[v]=1; c[v]=1;
            if (i%sshu[j]==0) {
                c[v]=c[i]+1;
                mu[v]=0; d[v]=d[i]*(c[v]+1)/c[v];
                break;
            }
            mu[v]=-mu[i]; d[v]=d[i]<<1;
        }
    }

    for (int i=1;i<N;i++)
        sm[i]=sm[i-1]+mu[i]*mu[i],
        sd[i]=sd[i-1]+d[i];
}

ll D(ll x) {
    if (x<N) return (ll)sd[x];
    ll ans=0,last;
    for (ll i=1;i<=x;i=last+1) {
        last=x/(x/i);
        ans+=(x/i*(ll)(last-i+1));
    }
    return ans;
}

ll Mu(ll x) {
    if (x<N) return sm[x];
    ll ans=0,last;
    for (ll i=1;i*i<=x;i++) 
        ans+=(mu[i]*(n/(i*i)));
    return ans;
}

ll solve(ll n) {
    ll ans=0,last;
    for (ll i=1;i<=n;i=last+1) { 
        last=n/(n/i);
        ans+=D(n/i)*(Mu(last)-Mu(i-1));
    }
    return ans;
}

int main()
{
    int T;
    prepare();
    scanf("%d",&T);
    while (T--) {
        scanf("%lld",&n);
        printf("%lld\n",solve(n));
    }
    return 0;
}