【BZOJ3529】数表(莫比乌斯反演,树状数组)

【BZOJ3529】数表(莫比乌斯反演,树状数组)

题解

首先不管\(A\)的范围的限制
要求的东西是
\[\sum_{i=1}^n\sum_{j=1}^m\sigma(gcd(i,j))\]
其中\(\sigma(x)\)表示\(x\)的约数之和

约数之和是一个积性函数,可以线性筛
具体的做法请参考皮皮亮的Blog

根据常见的套路
\(gcd\)给提出来
\[\sum_{d=1}^n\sigma(d)\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=d]\]
后面那个东西不用说了吧。。
我原来已经推过很多遍啦
可以参考原来推过的式子
所以我直接写啦
\[\sum_{d=1}^n\sigma(d)\sum_{i=1}^{n/d}\mu(i)[\frac{n}{id}][\frac{m}{id}]\]

还是一样的,设\(T=id\)
然后提出来
\[\sum_{T=1}^n[\frac{n}{T}][\frac{m}{T}]\sum_{d|T}\sigma(d)\mu(\frac{T}{d})\]

如果没有\(A\)的限制
这就可以直接筛后面的那个狄利克雷卷积
然后\(O(\sqrt n)\)的回答每一组询问

现在有了\(A\)的限制,貌似筛出来也没有什么用了
考虑到\(n,m<=10^5\)应该可以乱搞了
所以干脆就先筛出\(\sigma\)\(\mu\)

然后,这个限制还是很不好搞???
范围这么小
那我就暴力呀
离线询问后按照\(A\)排序
然后把\(\sigma\)排序
依次插到相对应的位置
因为要求前缀,就搞一个树状数组来算前缀和就行啦

最后是取膜,为啥要用\(long\ long\)
\(int\)自然溢出就行啦

暴力出奇迹

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
using namespace std;
#define MAX 100000
inline int read()
{
    int x=0,t=1;char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}
bool zs[MAX+1];
int pri[MAX],tot;
int mu[MAX+1],sig[MAX+1],sumd[MAX+1],powd[MAX+1];
int T,N,ans[MAX];
struct Ask{int n,m,A,id;}q[MAX];
int p[MAX+1],c[MAX+1];
bool operator<(Ask a,Ask b){return a.A<b.A;}
void pre(int N)
{
    zs[1]=true;sig[1]=mu[1]=p[1]=1;
    for(int i=2;i<=N;++i)
    {
        p[i]=i;
        if(!zs[i])pri[++tot]=i,mu[i]=-1,sig[i]=i+1,sumd[i]=i+1,powd[i]=i;
        for(int j=1;j<=tot&&i*pri[j]<=N;++j)
        {
            zs[i*pri[j]]=true;
            if(i%pri[j])
            {
                mu[i*pri[j]]=-mu[i];
                sig[i*pri[j]]=sig[i]*sig[pri[j]];
                sumd[i*pri[j]]=1+pri[j];
                powd[i*pri[j]]=pri[j];
            }
            else
            {
                powd[i*pri[j]]=powd[i]*pri[j];
                sumd[i*pri[j]]=sumd[i]+powd[i*pri[j]];
                sig[i*pri[j]]=sig[i]/sumd[i]*sumd[i*pri[j]];
                break;
            }
        }
    }
}
bool cmp(int a,int b){return sig[a]<sig[b];}
inline int lowbit(int x){return x&(-x);}
void Add(int x,int w){while(x<=N)c[x]+=w,x+=lowbit(x);}
int Query(int x){int ret=0;while(x)ret+=c[x],x-=lowbit(x);return ret;}
int Solve(int x,int y)
{
    int now=0,lst=0,ret=0;
    for(int i=1,j;i<=x;i=j+1)
    {
        j=min(x/(x/i),y/(y/i));
        now=Query(j);
        ret+=(x/i)*(y/i)*(now-lst);
        lst=now;
    }
    return ret;
}
int main()
{
    T=read();
    for(int i=1;i<=T;++i)
    {
        q[i].id=i,q[i].n=read(),q[i].m=read(),q[i].A=read();
        if(q[i].n>q[i].m)swap(q[i].n,q[i].m);
        N=max(N,q[i].n);
    }
    pre(N);
    sort(&q[1],&q[T+1]);sort(&p[1],&p[N+1],cmp);
    for(int i=1,j=1;i<=T;++i)
    {
        for(;j<=N&&sig[p[j]]<=q[i].A;++j)
            for(int k=p[j];k<=N;k+=p[j])
                if(mu[k/p[j]])Add(k,sig[p[j]]*mu[k/p[j]]);
        ans[q[i].id]=Solve(q[i].n,q[i].m);
    }
    for(int i=1;i<=T;++i)
    {
        if(ans[i]<0)ans[i]+=2147483647,ans[i]++;
        printf("%d\n",ans[i]);
    }
    return 0;
}

转载于:https://www.cnblogs.com/cjyyb/p/8303141.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值