bzoj3994:[SDOI2015]约数个数和

本文深入探讨了莫比乌斯反演在处理约数问题中的应用,详细解析了一个核心公式,展示了如何通过数学变换简化复杂求和过程,同时提供了高效算法实现,包括预处理和数论分块优化。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

传送门

莫比乌斯反演,但是约数有点难搞诶

有一个式子(想想挺显然的,可以保证每个约数只被算一次)
\[ d(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1] \]
然后就正常了对吧,设
\[ f(d)=\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{x|i}\sum_{y|j}[gcd(x,y)==d] \]

然后设
\[ g(n)=\sum_{n|d}f(d)=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{x|i}\sum_{y|j}[n|gcd(x,y)] \]
再化一下,可能跳跃的有点多,自行理解下
\[ g(n)=\sum_{i=1}^{N/n}\sum_{j=1}^{M/n}\lfloor \frac{N}{in} \rfloor\lfloor \frac{M}{jn} \rfloor \]
反演:
\[ f(d)=\sum_{d|n}\mu(\frac{n}{d})g(n)=\sum_{d|n}\mu(\frac{n}{d})\sum_{i=1}^{N/n}\sum_{j=1}^{M/n}\lfloor \frac{N}{in} \rfloor\lfloor \frac{M}{jn} \rfloor \]
\(T=\frac{n}{d}\)
\[ f(d)=\sum_{T=1}^{min(N,M)/d}\mu(T)\sum_{i=1}^{N/Td}\sum_{j=1}^{M/Td}\lfloor \frac{N}{iTd} \rfloor\lfloor \frac{M}{jTd} \rfloor \]
由于我们要求的\(ans=f(1)\)


\[ ans=f(1)=\sum_{T=1}^{min(N,M)}\mu(T)\sum_{i=1}^{N/T}\sum_{j=1}^{M/T}\lfloor \frac{N}{iT} \rfloor\lfloor \frac{M}{jT} \rfloor \]
如果直接数论分块的话,时间复杂度是\(O(TN)\)的,并不能通过

但是后面的\(\sum_{i=1}^{N/T}\sum_{j=1}^{M/T}\lfloor \frac{N}{iT} \rfloor\lfloor \frac{M}{jT} \rfloor\)可以考虑预处理,预处理时间复杂度\(O(n\sqrt{n})\)

总时间复杂度\(O(n\sqrt{n}+T\sqrt{n})\)

代码:

#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
void read(int &x) {
    char ch; bool ok;
    for(ok=0,ch=getchar(); !isdigit(ch); ch=getchar()) if(ch=='-') ok=1;
    for(x=0; isdigit(ch); x=x*10+ch-'0',ch=getchar()); if(ok) x=-x;
}
#define rg register
const int maxn=5e4+10;
int T,n,m,pri[maxn],tot,mu[maxn],sum[maxn];
bool vis[maxn];long long ans,f[maxn];
void prepare()
{
    mu[1]=1;
    for(rg int i=2;i<=5e4;i++)
    {
        if(!vis[i])pri[++tot]=i,mu[i]=-1;
        for(rg int j=1;j<=tot&&pri[j]*i<=5e4;j++)
        {
            vis[pri[j]*i]=1;
            if(!(i%pri[j]))break;
            else mu[i*pri[j]]=-mu[i];
        }
    }
    for(rg int i=1;i<=5e4;i++)sum[i]=sum[i-1]+mu[i];
}
void solve()
{
    for(rg int x=1;x<=5e4;x++)
    {
        long long ans=0;
        for(rg int i=1,j;i<=x;i=j+1)
        {
            j=x/(x/i);
            ans+=1ll*(x/i)*(j-i+1);
        }
        f[x]=ans;
    }
}
int main()
{
    read(T),prepare(),solve();
    while(T--)
    {
        read(n),read(m);if(n>m)swap(n,m);ans=0;
        for(rg int i=1,j;i<=n;i=j+1)
        {
            j=min(n/(n/i),m/(m/i));
            ans+=(sum[j]-sum[i-1])*f[n/i]*f[m/i];
        }
        printf("%lld\n",ans);
    }
}

转载于:https://www.cnblogs.com/lcxer/p/10544258.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值