[51Nod1584]加权约数和-线性筛-莫比乌斯反演

加权约数和

Description

去年的 tangjz 非常喜欢做数论题，但是一年以后的 tangjz 却不那么会做了。
在整理以前的试题时，他发现了这样一道题目：“求 ∑σ(i) ，其中 1≤i≤N ， σ(i) 表示 i 的约数之和。”
现在他长大了，题目也变难了，所以麻烦你来帮他解决一道数论题吧。
他需要你求如下表达式的值：
∑Ni=1∑Nj=1max(i,j)⋅σ(i⋅j)
其中 max(i,j) 表示 i 和 j 里的最大值， σ(i⋅j) 表示 i⋅j 的约数之和。
例如当 N=2 的时候，由 σ(1)=1,σ(2)=1+2=3,σ(4)=1+2+4=7 可知，答案应为 1⋅σ(1⋅1)+2⋅σ(1⋅2)+2⋅σ(2⋅1)+2⋅σ(2⋅2)=27 。
他发现答案有点大，所以你只需要告诉他答案模 1000000007 的值即可。

Input

每个测试点含有多组测试数据。
第一行是一个正整数T，表示接下来有T组测试数据。（1≤T≤50000）
接下来的T行，每组测试数据占一行。
每行有一个正整数N，含义如描述所示。（1≤N≤1000000）

Output

共有T行。对于每组测试数据，输出一行信息”Case #x: y”。
其中x表示对应的是第几组测试数据，y表示相应的答案模1000000007的值。

Input示例

5
1
2
3
4
5

Output示例

Case #1: 1
Case #2: 27
Case #3: 162
Case #4: 686
Case #5: 1741

---

看杜教筛结果被文末习题骗了过来……
推式子共用了10min多，实现方面调了2h+…….
这真的是一道推式子的题吗……

---

思路:
首先，考虑如何处理max的问题。
考虑想象一个矩阵，那么可以发现其max值沿对角线对称。
那么很自然的使用两倍的上三角矩阵减去对角线来代替原式:

$$ans=2\sum\limits_{i=1}^ni\sum\limits_{j=1}^i\sigma(i*j)-\sum\limits_{i=1}^ni*\sigma(i^2)$$

可以发现对角线部分可以线性筛预处理，那么只考虑前半部分，即

$$ans'=\sum\limits_{i=1}^ni\sum\limits_{j=1}^i\sigma(i*j)$$

考虑拆开$\sigma(i*j)$:

$$\sigma(i*j)=\sum\limits_{x|i}\sum\limits_{y|j}x*y*e(gcd(x,y))$$
反演:
$$\sigma(i*j)=\sum\limits_{x|i}\sum\limits_{y|j}x*y\sum\limits_{d|x,y}\mu(d)$$
将 $d$提前:
$$\sigma(i*j)=\sum\limits_{d|i,j}\mu(d)\sum\limits_{x|\frac{i}{d}}x\sum\limits_{y|\frac{i}{d}}y$$
将后面的两个求和符号化简:
$$\sigma(i*j)=\sum\limits_{d|i,j}\mu(d)\sigma(\frac{i}{d})\sigma(\frac{j}{d})$$

至此完成这一部分的化简!

于是带回原式:

$$ans'=\sum\limits_{i=1}^ni\sum\limits_{j=1}^i\sum\limits_{d|i,j}\mu(d)\sigma(\frac{i}{d})\sigma(\frac{j}{d})$$

将$d$提前:

$$ans'=\sum\limits_{d=1}^nd^2\mu(d)\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}i\sigma(i) \sum\limits_{j=1}^i \sigma(j)$$

带回总式子:

$$ans=2\sum\limits_{d=1}^nd^2\mu(d)\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}i\sigma(i) \sum\limits_{j=1}^i \sigma(j)-\sum\limits_{i=1}^ni*\sigma(i^2)$$

考虑$n$的范围不大，于是线性筛出如下内容:

$mu[n]=\mu(n)$
$d[n]=\sigma(n)$
$sumd[n]=\sum\limits_{i=1}^n \sigma(i)$
$sum[n]=\sum\limits_{i=1}^{n}i\sigma(i) \sum\limits_{j=1}^i \sigma(j)$ (代码中未出现)
$ds[n]=\sigma(n^2)$
$sumds[n]=\sum\limits_{i=1}^ni*\sigma(i^2)$
此处为全文精华部分，由于篇幅原因各种实现详见代码

于是现在就可以整除分块做到单词询问$O(\sqrt{n})$了——
然而这并没有什么卯月。询问数巨大，因此依旧会T。

考虑$O(n\ln(n))$预处理出$ans'$的答案数组$f[n]$，每次$O(1)$回答询问，这样就能过了~

这真的不是在考线性筛的各种操作吗……

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>

using namespace std;

typedef long long ll;
const int N=1e6+7;
const ll md=1e9+7;

ll n;
ll mu[N],pri[N>>1],ptop;
ll sumd[N],sumds[N];
ll mi1[N],mi2[N],d[N];
ll ds[N],mi1s[N],mi2s[N],f[N];
bool npri[N];

inline void init()
{
    mu[1]=ds[1]=d[1]=1;
    sumds[1]=sumd[1]=1;
    for(ll i=2;i<N;i++)
    {
        if(!npri[i])
        {
            pri[++ptop]=i;
            mu[i]=-1;

            d[i]=i+1;
            mi1[i]=i+1;
            mi2[i]=i;

            ds[i]=(ll)i*i+i+1ll;
            mi1s[i]=(ll)i*i+i+1ll;
            mi2s[i]=(ll)i*i;
        }
        for(int j=1;j<=ptop && i*pri[j]<N;j++)
        {
            npri[i*pri[j]]=1;
            if(i%pri[j])
            {
                mu[i*pri[j]]=-mu[i];

                d[i*pri[j]]=d[i]*d[pri[j]];
                mi1[i*pri[j]]=1+pri[j];
                mi2[i*pri[j]]=pri[j];

                ds[i*pri[j]]=ds[i]*ds[pri[j]];
                mi1s[i*pri[j]]=(ll)pri[j]*pri[j]+(ll)pri[j]+1ll;
                mi2s[i*pri[j]]=(ll)pri[j]*pri[j];
            }
            else
            {
                mu[i*pri[j]]=0;

                mi2[i*pri[j]]=mi2[i]*pri[j];
                mi1[i*pri[j]]=mi1[i]+mi2[i*pri[j]];
                d[i*pri[j]]=d[i]/mi1[i]*mi1[i*pri[j]];

                mi2s[i*pri[j]]=mi2s[i]*pri[j]*pri[j];
                mi1s[i*pri[j]]=mi1s[i]+mi2s[i]*((ll)pri[j]*pri[j]+pri[j]);
                ds[i*pri[j]]=ds[i]/mi1s[i]*mi1s[i*pri[j]];
                break;
            }
        }
        sumd[i]=(sumd[i-1]+d[i])%md;
        sumds[i]=(sumds[i-1]+ds[i]*i%md)%md;
    }

    for(int i=1;i<N;i++)
        if(mu[i])
        {
            for(int j=i;j<N;j+=i)
            {
                if(mu[i]==1)
                    f[j]+=(ll)i*j%md*sumd[j/i]%md*d[j/i]%md;
                else
                    f[j]+=md-(ll)i*j%md*sumd[j/i]%md*d[j/i]%md;
                if(f[j]>=md)
                    f[j]-=md;
            }
        }
    for(int i=1;i<N;i++)
    {
        f[i]+=f[i-1];
        if(f[i]>=md)
            f[i]-=md;
    }
}

int main()
{
    init();
    int t;
    scanf("%d",&t);
    for(int tim=1;tim<=t;tim++)
    {
        scanf("%lld",&n);
        printf("Case #%d: %lld\n",tim,((2ll*f[n]%md-sumds[n])%md+md)%md);
    }

    return 0;
}