51nod1584 加权约数和

我们把答案写成$2\sum_{i=1}^n\sum_{j=1}^ii\sigma(ij)-\sum_{i=1}^n\sigma(i^2)$，对两部分分别求和。
后一部分比较简单，$\sigma(i^2)$是积性函数，且$\sigma((p^xq\cdot p)^2)=(\sigma((p^x)^2)+p^{2x+1}+p^{2x+2})\cdot \sigma(q^2)$，可以线性筛。
考虑前一部分，首先我们知道，$\sigma(mn)=\sum_{i|n}\sum_{j|m}ij e(\gcd(\frac n i,j))$。那么

$$\begin{align} & \sum_{i=1}^n i\sum_{j=1}^i\sigma(ij) \\ = & \sum_{i=1}^ni\sum_{j=1}^i\sum_{a|i}a\sum_{b|j}\frac j b\sum_{k|\gcd(a,b)} \mu(k) \\ = & \sum_{k=1}^n\mu(k)\sum_{i=1}^{\lfloor\frac n k\rfloor}ik\sum_{a|i}ak\sum_{j=1}^i \sum_{b|j}\frac j b \\ = & \sum_{k=1}^n\mu(k)k^2\sum_{i=1}^{\lfloor\frac n k\rfloor} g(i) \\ = & \sum_{K=1}^n\sum_{k|K}\mu(k)k^2g(\frac K k) \\ \end{align}$$
在求出 $g$函数之后，我们就可以$O(n\log n)$计算出 $1$到$n$的答案。而 $g(n) =n\sigma(n)\sum_{i=1}^n\sigma(i)$。线性筛求出 $\sigma$的前缀和就可以了。
最后的复杂度是 $O(n\log n+q)$。

#include<cstdio>
#include<algorithm>
using namespace std;
#define LL long long
const int p=1000000007,maxn=1000000;
int inv[maxn+10],prm[maxn+10],mu[maxn+10],mnp[maxn+10],mnpw[maxn+10],
f[maxn+10],g[maxn+10],h[maxn+10],
ssigma2[maxn+10],sigma[maxn+10],
n,tot;
int inc(int x,int y)
{
    x+=y;
    return x>=p?x-p:x;
}
int dec(int x,int y)
{
    x-=y;
    return x<0?x+p:x;
}
int pow(int b,int k)
{
    int ret=1;
    for (;k;k>>=1,b=(LL)b*b%p)
        if (k&1) ret=(LL)ret*b%p;
    return ret;
}
void init()
{
    int x;
    inv[1]=mu[1]=sigma[1]=ssigma2[1]=1;
    for (int i=2;i<=maxn;i++)
    {
        if (!mnp[i])
        {
            prm[++tot]=i;
            mnp[i]=mnpw[i]=i;
            mu[i]=p-1;
            inv[i]=pow(i,p-2);
            sigma[i]=i+1;
            ssigma2[i]=inc(inc((LL)i*i%p,i),1);
        }
        for (int j=1;j<=tot&&(LL)i*prm[j]<=maxn;j++)
        {
            x=i*prm[j];
            inv[x]=(LL)inv[i]*inv[prm[j]]%p;
            if (i%prm[j])
            {
                mnp[x]=mnpw[x]=prm[j];
                mu[x]=dec(0,mu[i]);
                sigma[x]=(LL)sigma[i]*sigma[prm[j]]%p;
                ssigma2[x]=(LL)ssigma2[i]*ssigma2[prm[j]]%p;
            }
            else
            {
                mnp[x]=prm[j];
                mnpw[x]=mnpw[i]*prm[j];
                mu[x]=0;
                sigma[x]=inc(sigma[i],(LL)prm[j]*mnpw[i]%p*sigma[i/mnpw[i]]%p);
                ssigma2[x]=inc(ssigma2[i],(LL)inc((LL)mnpw[i]*mnpw[i]%p*mnp[i]%p,(LL)mnpw[i]*mnpw[i]%p*mnp[i]%p*mnp[i]%p)*ssigma2[i/mnpw[i]]%p);
                break;
            }
        }
    }
    for (int i=1;i<=maxn;i++)
    {
        ssigma2[i]=inc(ssigma2[i-1],(LL)i*ssigma2[i]%p);
        h[i]=inc(h[i-1],sigma[i]);
        g[i]=(LL)i*h[i]%p*sigma[i]%p;
        for (int j=i;j<=maxn;j+=i) f[j]=inc(f[j],(LL)g[i]*mu[j/i]%p*(j/i)%p*(j/i)%p);
        f[i]=inc(f[i-1],f[i]);
    }
}
int main()
{
    init();
    int T;
    scanf("%d",&T);
    for (int K=1;K<=T;K++)
    {
        scanf("%d",&n);
        printf("Case #%d: %d\n",K,dec(inc(f[n],f[n]),ssigma2[n]));
    }
}