【51Nod1584】加權約數和（數論）

【51Nod1584】加權約數和（數論）

題面

51Nod

題解

要求的是\[\sum_{i=1}^n\sum_{j=1}^n max(i,j)\sigma(ij)\]
這個\(max\)太討厭了，直接枚舉一半乘個二。
\[2\sum_{i=1}^n\sum_{j=1}^{i}i\sigma(ij)-\sum_{i=1}^ni\sigma(i^2)\]
後面這一半能夠直接預處理，只須要把\(i\)分解，能夠作到調和級數的複雜度。
只考慮前面這一半，顯然只須要考慮的是\(\sigma(ij)\)這個東西。
那麼咱們考慮在\(i\)中枚舉一個約數，在\(j\)中枚舉一個約數，而後把這兩個約數合併一下，看看能不能讓每一個約數只被計算一次。
\[\sigma(ij)=\sum_{u|i}\sum_{v|j}[gcd(u,\frac{j}{v})=1]uv\]
證實的話，大概就是咱們的目標是讓每一個約數只被計算一次，首先在\(i\)中枚舉一個約數確定沒有問題，在\(j\)中枚舉一個質因數也沒有問題。對於\(uv\)這個數而言，咱們把只在\(i\)中有的因子和只在\(j\)中有的因子給丟掉，只考慮在\(i,j\)中都含有的因子\(u',v'\)，對於一個數\(uv\)而言，可能算重的狀況是\(u'\)從\(v'\)那裏搶走了一個質因子，而此時\(\frac{j}{v}\)就會對應的乘上那個質因子，使得\(gcd\neq 1\)，因此每一個數只會被計算一次。
有了這個式子就很好搞了，首先把這個式子換一個形式：
\[\sigma(ij)=\sum_{u|i}\sum_{v|j}[gcd(u,v)=1]\frac{uj}{v}\]
帶回去獲得：
\[\sum_{i=1}^n i\sum_{j=1}^i\sum_{u|i}\sum_{v|j}[gcd(u,v)=1]\frac{uj}{v}\]
考慮對於每個\(i\)分別計算答案，因此咱們設
\[\begin{aligned} f[n]&=n\sum_{j=1}^n\sum_{u|n}\sum_{v|j}[gcd(u,v)=1]\frac{uj}{v}\\ &=n\sum_{j=1}^n \sum_{u|n}\sum_{v|j}\frac{uj}{v}\sum_{k|u,k|v}\mu(k)\\ &=n\sum_{j=1}^n\sum_{k|n,k|j}\mu(k)\sum_{k|u,u|n}\sum_{k|v,v|j}\frac{uj}{v}\\ &=n\sum_{j=1}^n\sum_{k|n,k|j}\mu(k)(k\sigma(\frac{n}{k}))\sigma(\frac{j}{k})\\ &=n\sum_{k|n}\mu(k)k\sigma(\frac{n}{k})\sum_{i=1}^{n/k}\sigma(i) \end{aligned}\]
而後就是前綴和計算就好了。
全部東西能夠線性篩，中間要求逆就直接快速冪了。。。。

#include<iostream>
#include<cstdio>
using namespace std;
#define ll long long
#define MOD 1000000007
#define MAX 1000100
inline int read()
{
    int x=0;bool t=false;char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=true,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return t?-x:x;
}
int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
bool zs[MAX];
int pri[MAX],tot;
int mu[MAX],sig[MAX],ssig[MAX],pw[MAX],spw[MAX],dpw[MAX],dspw[MAX],dsig[MAX];
int f[MAX],s[MAX];
void Sieve(int n)
{
    mu[1]=1;dsig[1]=sig[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!zs[i])
        {
            pri[++tot]=i,mu[i]=MOD-1;
            sig[i]=i+1,pw[i]=i,spw[i]=i+1;
            dsig[i]=dspw[i]=(1+i+1ll*i*i)%MOD;dpw[i]=1ll*i*i%MOD;
        }
        for(int j=1;j<=tot&&i*pri[j]<=n;++j)
        {
            zs[i*pri[j]]=true;
            if(i%pri[j])
            {
                mu[i*pri[j]]=MOD-mu[i];
                sig[i*pri[j]]=1ll*sig[i]*sig[pri[j]]%MOD;
                pw[i*pri[j]]=pri[j],spw[i*pri[j]]=1+pri[j];
                dsig[i*pri[j]]=1ll*dsig[i]*dsig[pri[j]]%MOD;
                dpw[i*pri[j]]=dpw[pri[j]];
                dspw[i*pri[j]]=dspw[pri[j]];
            }
            else
            {
                mu[i*pri[j]]=0;
                sig[i*pri[j]]=1ll*sig[i]*fpow(spw[i],MOD-2)%MOD*(spw[i]+pw[i]*pri[j])%MOD;
                pw[i*pri[j]]=pw[i]*pri[j];
                spw[i*pri[j]]=(spw[i]+pw[i]*pri[j])%MOD;
                dspw[i*pri[j]]=(dspw[i]+1ll*dpw[i]*pri[j]%MOD+1ll*dpw[i]*pri[j]%MOD*pri[j]%MOD)%MOD;
                dsig[i*pri[j]]=1ll*dsig[i]*fpow(dspw[i],MOD-2)%MOD*dspw[i*pri[j]]%MOD;
                dpw[i*pri[j]]=1ll*dpw[i]*pri[j]%MOD*pri[j]%MOD;
                break;
            }
        }
    }
    for(int i=1;i<=n;++i)dsig[i]=1ll*dsig[i]%MOD*i%MOD;
    for(int i=1;i<=n;++i)ssig[i]=(ssig[i-1]+sig[i])%MOD;
    for(int i=1;i<=n;++i)
    {
        if(mu[i])
            for(int j=i;j<=n;j+=i)
                f[j]=(f[j]+1ll*mu[i]*i%MOD*sig[j/i]%MOD*ssig[j/i])%MOD;
        f[i]=1ll*f[i]*i%MOD;s[i]=(s[i-1]+2ll*f[i]+MOD-dsig[i])%MOD;
    }
}
int main()
{
    Sieve(MAX-1);
    int T=read();
    for(int i=1;i<=T;++i)
        printf("Case #%d: %d\n",i,s[read()]);
    return 0;
}