51nod 1584 加权约数和 莫比乌斯反演

题意

去年的 tangjz 非常喜欢做数论题，但是一年以后的 tangjz 却不那么会做了。
在整理以前的试题时，他发现了这样一道题目：“求 ∑σ(i) ，其中 1≤i≤N ， σ(i) 表示 i 的约数之和。”
现在他长大了，题目也变难了，所以麻烦你来帮他解决一道数论题吧。
他需要你求如下表达式的值：
$\sum_{i=1}^n\sum_{j=1}^nmax(i,j)*\sigma(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 的值即可。
1≤T≤50000,1≤N≤1000000

分析

考虑先把那个讨厌的max去掉。

$$ans=2*\sum_{i=1}^ni\sum_{j=1}^i\sigma(ij)-\sum_{i=1}^ni*\sigma(i^2)$$
来继续推式子
$$sum=\sum_{i=1}^ni\sum_{j=1}^i\sigma(ij)$$
$$=\sum_{i=1}^ni\sum_{j=1}^i\sum_{p|i}\sum_{q|j}[(p,q)=1]\frac{pj}{q}$$
反演一下可以得到 $$sum=\sum_{d=1}^n\mu(d)d^2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i\sum_{j=1}^i\sum_{p|i}\sum_{q|j}pq$$
$$=\sum_{d=1}^n\mu(d)d^2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i*\sigma(i)\sum_{j=1}^n\sigma(j)$$
这已经是一个分块的形式了，但在这题仍然会T。我们考虑继续推式子。
设 $T=id$可以得到 $$sum=\sum_{T=1}^n\sum_{d|T}\mu(d)*d*T*\sigma(\frac{T}{d})*\sum_{i=1}^{\frac{T}{d}}\sigma(i)$$
然后就可以 $O(nln(n))$预处理答案 $O(1)$询问了。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;

const int N=1000005;
const int MOD=1000000007;

int sig[N],mu[N],s[N],sig2[N],tot,prime[N],low[N],tmp[N],ans[N];
bool not_prime[N];

int read()
{
    int x=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}

void get_prime(int n)
{
    sig[1]=mu[1]=sig2[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++tot]=i,sig[i]=i+1,mu[i]=-1,low[i]=i,sig2[i]=(1+i+(LL)i*i%MOD)%MOD;
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                low[i*prime[j]]=low[i]*prime[j];
                sig[i*prime[j]]=((LL)sig[i]*prime[j]+(LL)sig[i/low[i]])%MOD;
                sig2[i*prime[j]]=((LL)sig2[i]*prime[j]%MOD*prime[j]+(LL)sig2[i/low[i]]+(LL)sig2[i/low[i]]*prime[j])%MOD;
                break;
            }
            low[i*prime[j]]=prime[j];
            sig[i*prime[j]]=(LL)sig[i]*sig[prime[j]]%MOD;
            sig2[i*prime[j]]=(LL)sig2[i]*sig2[prime[j]]%MOD;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for (int i=1;i<=n;i++) mu[i]+=mu[i]<0?MOD:0,s[i]=(s[i-1]+sig[i])%MOD,tmp[i]=((LL)tmp[i-1]+(LL)i*sig2[i])%MOD;
    for (int i=1;i<=n;i++)
        for (int j=i;j<=n;j+=i)
        {
            ans[j]+=(LL)mu[i]*i%MOD*j%MOD*sig[j/i]%MOD*s[j/i]%MOD;
            ans[j]-=ans[j]>=MOD?MOD:0;
        }
    for (int i=1;i<=n;i++) ans[i]+=ans[i-1],ans[i]-=ans[i]>=MOD?MOD:0;
}

int main()
{
    get_prime(1000000);
    int T=read();
    for (int i=1;i<=T;i++)
    {
        int x=read();
        printf("Case #%d: %d\n",i,(ans[x]*2%MOD+MOD-tmp[x])%MOD);
    }
    return 0;
}