bzoj3529: [Sdoi2014]数表

链接

　　http://www.lydsy.com/JudgeOnline/problem.php?id=3529

题解

　　令$g(i)$表示$i$的所有因数之和，这个可以$O(Nlog_2N)$算出来。
　　$f(i)$表示$x\in[1,N],y\in[1,M]$中有多少数对$(x,y)$满足$gcd(x,y)=i$。如果不考虑$\le a$的限制，那就成了求表格里所有数的和，那就要枚举最大公约数$x$，看有多少数对的$gcd$等于$x$，数对的个数乘上$x$，然后所有的相加即可得答案，即
　　

$$ans=\sum_x^ng(x)f(x)$$
反演
　　 $$ans=\sum_{x=1}^ng(x)\sum_{x|d}^n\mu(\frac{d}{x})\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$$
　　这样是没法优化的，这里有一个技巧，我把它叫做 改变枚举顺序，即有原来的先枚举 $x$再枚举其倍数$d$，变为先枚举 $d$，再枚举$d$的因子 $x$。
　　 $$ans=\sum_{d=1}^n\sum_{x|d}g(x)\mu(\frac{d}{x})\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$$
　　当 $d$确定时，由于我们没有考虑$\le a$的限制，所以 $x$取那些值也确定了。因此，对于每一个$d$，预处理出 $s[d]=\sum_{x|d}^ng(x)\mu(\frac{d}{x})$，式子就成了
　　 $$ans=\sum_{d=1}^ns[d]\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$$
　　就可以分块统计答案了，时间复杂度 $O(T\sqrt N)$
　　可是…有 $\le a$这个限制。可以把 $g(x)>0$的数的 $g(x)$暂时变成0，暴力重构前缀和，这样对时间的牺牲很大。
　　离线，把所有询问按照 $a$排个序，所有的数按照$g(x)$排个序，搞两个指针，每次把合法的 $x$加入，用树状数组维护前缀和。
　　最终时间复杂度$O(Nlog_2N+T\sqrt{N}+Tlog_2T)$

代码

//莫比乌斯反演
#include <cstdio>
#include <algorithm>
#define maxn 100000
#define ll long long
#define lowbit(x) (x&-x)
#define mod ((ll)1<<31)
using namespace std;
ll prime[maxn+10], mark[maxn+10], mu[maxn+10], s[maxn+10], g[maxn+10], numg[maxn+10],
   quiry[maxn][4], T, numq[maxn];
void add(ll pos, ll v){for(;pos<=maxn;pos+=lowbit(pos))s[pos]+=v;}
ll sum(ll pos)
{
    ll ans=0;
    for(;pos;pos-=lowbit(pos))ans+=s[pos];
    return ans;
}
bool cmp1(ll a, ll b){return g[a]<g[b];}
bool cmp2(ll a, ll b){return quiry[a][3]<quiry[b][3];}
void init()
{
    ll i, j;
    mu[1]=1;
    for(i=2;i<=maxn;i++)
    {
        if(!mark[i])prime[++prime[0]]=i,mu[i]=-1;
        for(j=1;j<=prime[prime[0]] and i*prime[j]<=maxn;j++)
        {
            mark[i*prime[j]]=1;
            if(i%prime[j]==0){mu[i*prime[j]]=0;break;}
            mu[i*prime[j]]=-mu[i];
        }
    }
    for(i=1;i<=maxn;i++)
        for(j=i;j<=maxn;j+=i)g[j]+=i;
    for(i=1;i<=maxn;i++)numg[i]=i;
    sort(numg+1,numg+maxn+1,cmp1);
    scanf("%lld",&T);
    for(i=1;i<=T;i++)scanf("%lld%lld%lld",quiry[i]+1,quiry[i]+2,quiry[i]+3);
    for(i=1;i<=T;i++)numq[i]=i;
    sort(numq+1,numq+T+1,cmp2);
}
ll calc(ll n, ll m)
{
    ll i, last, ans=0;
    if(n>m)swap(n,m);
    for(i=1;i<=n;i=last+1)
    {
        last=min(n/(n/i),m/(m/i));
        ans+=(sum(last)-sum(i-1))*(n/i)*(m/i);
    }
    return ans;
}
void work()
{
    ll i, p=1, j;
    for(i=1;i<=T;i++)
    {
        for(;g[numg[p]]<=quiry[numq[i]][3] and p<=maxn;p++)
            for(j=numg[p];j<=maxn;j+=numg[p])add(j,g[numg[p]]*mu[j/numg[p]]);
        quiry[numq[i]][0]=calc(quiry[numq[i]][1],quiry[numq[i]][2]);
    }
}
int main()
{
    init();
    work();
    for(ll i=1;i<=T;i++)printf("%lld\n",quiry[i][0]%mod);
    return 0;
}