51nod1584：加权约数和 （莫比乌斯反演）

题目传送门：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1584

---

题目分析：这题应该算是[SDOI2015]约数个数和的加强版吧，如果还没有做过这道题的同学可以先想一想这题。至于它的解法你们可以自行百度（或看看我写的题解）
好吧，回到这题。
我们要搞清楚一个东西：$\sigma_1(ij)$（即i*j的约数和）究竟等于什么？（以下部分是我拿到这题时自己YY出来的）

---

在[SDOI2015]约数个数和这题中，我们有一个很巧妙的公式：$\sigma_0(ij)=\sum_{x|i}\sum_{y|j}[(x,y)=1]$，它主要是由$k+q+1=\sum_{i=0}^{k}\sum_{j=0}^{q}[(p^i,p^j)=1]$推来的。那我们来看看它的本质是什么，我们做出$i=1->n,j=1->n$时$(p^i,p^j)$的表格：

$(p^i,p^j)$	$i=0$	$i=1$	……	$i=k$
$j=q$	$1$	$p$	……	$p^{min(k,q)}$
……	……	……	……	……
$j=1$	$1$	$p$	……	$p$
$j=0$	$1$	$1$	……	$1$

很明显，我们发现表格中1的个数刚好就是$k+q+1$，而表格中1的要求是i或j至少有一个是0。然而这种限制是不可能用于多个质数$p_i$合并的，所以我们才写成了$\sum_{i=0}^{k}\sum_{j=0}^{q}[(p^i,p^j)=1]$这种形式，这样所有的质数合并到一起时，就会变成$\sum_{x|i}\sum_{y|j}[(x,y)=1]$这种方便处理的形式。
现在我们要求的是$\sigma_1(ab)$。如果单独考虑一个质数$p$，令$a=p^k*……,b=p^q*……$，那么我们要在答案中贡献$(1+p+p^2+……+p^{k+q})$这个因子。由于列出gcd和lcm的表格，格子中都不可能出现$p^{k+q}$，于是我列出了$p^i*p^j$的表格：

$p^i*p^j$	$i=0$	$i=1$	……	$i=k$
$j=q$	$p^q$	$p^{q+1}$	……	$p^{k+q}$
……	……	……	……	……
$j=1$	$p$	$p^2$	……	$p^{k+1}$
$j=0$	$1$	$p$	……	$p^k$

这样我们就看得很清楚了，我们要求的是$\sum_{i=0}^{k}\sum_{j=0}^{q}p^{i+j}[i=0或j=q]$，也就是表格倒”L”字形的部分。然而这样还是没办法用于多个质数合并，于是我改了改，写成$\sum_{i=0}^{k}\sum_{j=0}^{q}p^{i+j}[(p^i,\frac{p^q}{p^j})=1]$，然后多个质数合并在一起之后，就变成了最终式子：

$$\sigma_1(ij)=\sum_{x|i}\sum_{y|j}xy[(x,\frac{j}{y})=1]$$

---

以上就是我在草稿纸上自己YY的过程，后来我才知道就是这样做的。
然后我们把这个式子带进原题中，基本上就求解了。
说实话，在我写这篇总结之前，网上已经有大神给出了十分优秀的题解，但鉴于这提示本人自己推导的，我还是将自己的推导过程写在下面，以作纪念吧。

---

题目要求：

$$\sum_{i=1}^{n}\sum_{j=1}^{n}max(i,j)\sigma_1(ij)$$
我们发现这个max很讨厌，于是令j只枚举到i，变形成：
$$2\sum_{i=1}^{n}\sum_{j=1}^{i}i\sigma_1(ij)-\sum_{i=1}^{n}i\sigma_1({i^2})$$
其中右边是要减去重复计算的ans(i,i)，这一部分用线性筛后再做一次前缀和即可O(1)得到。线性筛时若i%prime[j]!=0，直接用积性函数的性质$\sigma_1(i*prime[j])=\sigma_1(i)*\sigma_1(prime[j])$即可。若i是prime[j]的倍数，则$\sigma_1(i*prime[j])=\sigma_1(i)+\sigma_1(\frac{i}{prime[j]^{h(i)}})*(prime[j]^{h(i)+2}+prime[j]^{h(i)+1})$，其中h(i)是i的最小质因数的最高次幂，这个也可以用线性筛维护（唉，老套路了）。

---

然后考虑左边2后面那一坨，转化成：

$$\sum_{i=1}^{n}i\sum_{j=1}^{i}\sum_{x|i}\sum_{y|j}x\frac{j}{y}[(x,y)=1]$$
将x,y的枚举拉到i,j前面，重新令$a=\frac{i}{x},b=\frac{j}{y}$，整理得到：
$$\sum_{x=1}^{n}x^2\sum_{a=1}^{\lfloor\frac{n}{x}\rfloor}a\sum_{y=1}^{ax}\sum_{b=1}^{\lfloor\frac{ax}{y}\rfloor}b\sum_{k|x,k|y}\mu(k)$$
然后我们将k的枚举拉到最前面，重新令$i=\frac{x}{k},j=\frac{y}{k}$，得到：
$$\sum_{k=1}^{n}k^2\mu(k)\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}i^2\sum_{a=1}^{\lfloor\frac{n}{ik}\rfloor}a\sum_{y=1}^{ai}\sum_{b=1}^{\lfloor\frac{ai}{y}\rfloor}b$$
我们记$G(n)=\sum_{i=1}^{n}\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}d$（很明显这就是约数和的前缀和，可以用线性筛求），$F(n)=nG(n)\sum_{d|n}d$（很明显，这也是可以线性求的），这样原式子就变成了：
$$\sum_{k=1}^{n}k^2\mu(k)\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}F(i)$$
对F函数做前缀和，就可以做到预处理O(n)，询问O($\sqrt{n}$)，然而50000组数据，还是会超。

---

式子推到上面我就推不下去了，感觉没办法让时间更优。
后来上网看了看别人的题解，原来再换个枚举顺序就可以了……（我是智障）
记d=i*k，外层枚举d，得到：

$$\sum_{d=1}^{n}\sum_{k|d}k^2\mu(k)F(\frac{d}{k})$$
然后就变成了预处理O(n*ln(n))，询问O(1)。

---

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=1000010;
const long long M=1000000007;
typedef long long LL;

int pk[maxn];
LL h[maxn];
LL p[maxn];
LL mu[maxn];

bool vis[maxn];
int prime[maxn];
int cur=0;

LL g[maxn];
LL f[maxn];

int t,n;

void Preparation()
{
    h[1]=p[1]=mu[1]=1;
    for (int i=2; i<maxn; i++)
    {
        if (!vis[i]) pk[i]=i,h[i]=i+1,p[i]=i,p[i]*=i,p[i]=(p[i]+i+1)%M,mu[i]=-1,prime[++cur]=i;
        for (int j=1; j<=cur && i*prime[j]<maxn; j++)
        {
            int k=i*prime[j];
            vis[k]=true;
            if (i%prime[j]) pk[k]=prime[j],h[k]=h[i]*h[ pk[k] ],p[k]=p[i]*p[ pk[k] ]%M,mu[k]=-mu[i];
            else
            {
                pk[k]=pk[i]*prime[j];
                h[k]=h[i]+h[ i/pk[i] ]*pk[k];
                p[k]=pk[i]+pk[k],p[k]=p[k]*pk[k]%M*p[ i/pk[i] ]%M,p[k]=(p[k]+p[i])%M;
                mu[k]=0;
                break;
            }
        }
    }
    for (int i=1; i<maxn; i++) p[i]=p[i]*i%M;
    for (int i=2; i<maxn; i++) p[i]=(p[i]+p[i-1])%M;
    for (int i=1; i<maxn; i++) mu[i]=(mu[i]+M)*i%M*i%M;
    g[1]=1;
    for (int i=2; i<maxn; i++) g[i]=(g[i-1]+h[i])%M;
    for (int i=1; i<maxn; i++) g[i]=g[i]*i%M*h[i]%M;
    for (int i=1; i<maxn; i++)
        for (int j=1; i*j<maxn; j++)
            f[i*j]=(f[i*j]+ mu[i]*g[j]%M )%M;
    for (int i=2; i<maxn; i++) f[i]=(f[i]+f[i-1])%M;
}

int main()
{
    freopen("c.in","r",stdin);
    freopen("c.out","w",stdout);

    Preparation();
    scanf("%d",&t);
    for (int i=1; i<=t; i++)
    {
        scanf("%d",&n);
        LL ans=(f[n]*2LL-p[n]+M)%M;
        printf("Case #%d: %lld\n",i,ans);
    }

    return 0;
}