【BZOJ3994】约数个数和(SDOI2015)-莫比乌斯反演+数论分块

测试地址:约数个数和
做法:本题需要用到莫比乌斯反演+数论分块。
首先,因为ij的每个因数都可以唯一表示成pjq的形式,其中p|i,q|j,gcd(p,q)=1,所以有:
d(ij)=x|iy|j[gcd(x,y)=1]
把这个式子带进原式中去,有:
ans=i=1nj=1mx|iy|j[gcd(x,y)=1]
x,y换出来,得到:
ans=x=1ny=1m[gcd(x,y)=1]nxmy
由莫比乌斯反演定理,有d|nμ(d)=[n=1](可以从狄利克雷卷积的形式来理解,因为I=Ie,所以有e=μI),所以:
ans=x=1nnxy=1mmyd|gcd(x,y)μ(d)
d换出来,得到:
ans=d=1min(n,m)μ(d)x=1ndnxdy=1mdmyd
我们知道nxy=nxy,简单证明如下:
nxy=k,则有n=k(xy)+r,其中0r<xy,又令r=r0x+r1,其中0r1<x,可以知道0r0<y,那么有nx=ky+r0,所以nxy=k
所以上式可以写成:
ans=d=1min(n,m)μ(d)x=1ndndxy=1mdmdy
那么我们只需预处理出μ的前缀和,以及找到计算后面那一串东西的快速方法,就可以数论分块了。令S(n)=i=1nni,那么上式可以写成:
ans=d=1min(n,m)μ(d)S(nd)S(md)
我们又知道S(n)=i=1nd(i),为什么呢?因为看S的第一种表示,实际上是枚举一个数,这个数对答案做出它的倍数个数的贡献,那么就可以换乘枚举每个数,这个数对答案做出它的因数个数的贡献,所以两种表示是等价的。那么显然我们可以线性筛出d来预处理S,那么这个问题就解决了。
于是采用数论分块算ans,时间复杂度为O(Tn)
以下是本人代码:

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int T;
ll n[50010],m[50010],maxn=0;
ll prime[50010],mu[50010],summu[50010];
ll id[50010],d[50010],sumd[50010];
bool vis[50010]={0};

void calc(ll limit)
{
    mu[1]=d[1]=1;
    prime[0]=0;
    for(ll i=2;i<=limit;i++)
    {
        if (!vis[i])
        {
            prime[++prime[0]]=i;
            mu[i]=-1;
            id[i]=2;
            d[i]=2;
        }
        for(ll j=1;j<=prime[0]&&i*prime[j]<=limit;j++)
        {
            vis[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                id[i*prime[j]]=id[i]+1;
                d[i*prime[j]]=d[i]*(id[i]+1ll)/id[i];
                break;
            }
            mu[i*prime[j]]=-mu[i];
            id[i*prime[j]]=2;
            d[i*prime[j]]=d[i]*id[i*prime[j]];
        }
    }
    summu[0]=sumd[0]=0;
    for(ll i=1;i<=limit;i++)
    {
        summu[i]=summu[i-1]+mu[i];
        sumd[i]=sumd[i-1]+d[i];
    }
}

ll solve(ll n,ll m)
{
    ll ans=0;
    for(ll i=min(n,m);i>=1;i=max(n/(n/i+1),m/(m/i+1)))
    {
        ll l=max(n/(n/i+1),m/(m/i+1))+1,r=i;
        ans+=(summu[r]-summu[l-1])*sumd[n/i]*sumd[m/i];
    }
    return ans;
}

int main()
{
    scanf("%d",&T);
    for(int i=1;i<=T;i++)
    {
        scanf("%lld%lld",&n[i],&m[i]);
        maxn=max(maxn,max(n[i],m[i]));
    }

    calc(maxn);
    for(int i=1;i<=T;i++)
        printf("%lld
",solve(n[i],m[i]));

    return 0;
}
原文地址:https://www.cnblogs.com/Maxwei-wzj/p/9793408.html