【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;
}
posted @ 2019-07-06 09:51 小蒟蒻yyb 阅读(...) 评论(...) 编辑 收藏