BZOJ2301 Problem B 【莫比乌斯反演】【分块】

题目链接:https://www.lydsy.com/JudgeOnline/problem.php?id=2301

题解:
对于 axbcyd,这个条件,我们发现比较难以处理,这时候我们可以利用二维前缀和的思想,记 xbyd 时的答案为 A[b][d],那么Ans=A[b][d]A[b][c1]A[a1][d]+A[a1][c1],这样就把问题简化成了两个条件的问题。
下面推一下式子:
Ans=Σi=1b Σj=1d (gcd(i,j)==k)

=Σi=1b Σj=1d (gcd(i/k,j/k)==1)

我们发现 gcd(X,Y)==1,可以用卷积的单位元来做。

=Σi=1b Σj=1d (e(gcd(i/k,j/k))==1)

=Σi=1b Σj=1d ΣT|(i/k),T|(j/k)μ(T)

=ΣT|(i/k),T|(j/k) μ(T) Σi=1b/(kT) Σj=1d/(kT)1

=ΣT|(i/k),T|(j,k) μ(T)(b/(kT))(d/(kT))

对于(b/(kT))(d/(kT)),对于(b/k)(d/k),将其看做nm,发现可以分块即可。

原因:(b/(kT))=((b/k)/T) –> 将 b/k 看做 n,一项即为 n/T
可以从小到大枚举 T,相同的块可以直接跳过。
现在我们的任务变成了:如何求一个块的最大元素。这个严格证明参考了某大神的blog。
设这个块的 Ti
n=pi+q,那么对于 i,若 i++,则 q 会相应减少 p,而 q 最多减少 q/pp
所以此时的 i=i+q/p=(pi+q)/p,即为 n/(n/i),每次分块做即可
PS:小心爆 int,如果不分块的话极限数据是50000^2级别的,会TLE。
关于求μ()函数,由于其为积性函数,所以可以用线性筛解决,详见代码。

代码:

// by Balloons
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define mpr make_pair
#define debug() puts("okkkkkkkk")
#define rep(i,a,b) for(int (i)=(a);(i)<=(b);(i)++)

using namespace std;

typedef long long LL;

const int inf = 1 << 30;

const int maxn=50000;
int notpm[maxn+5],pm[maxn+5],mu[maxn+5],pcnt=0;
int tes;
int a,b,c,d,k;

/*
推导之后,我们发现就是一个式子:
sigma T=1 to min(x/k,y/k) miu(T) * (b/(k*T)) * (d/(k*T))
对于后面两项,我们发现可以分块。
原因:(b/(k*T))=((b/k)/T) --> 将 b/k 看做 n,一项即为 n/T,
可以从小到大枚举 T,相同的块可以直接跳过。
现在我们的任务变成了:如何求一个块的最大元素。设这个块的 T为i 
设 n=pi+q,那么对于 i,若 i++,则 q 会相应减少 p,而 q 最多减少 q/p个p,
所以此时的 i'=i+q/p=(pi+q)/p,即为 n/(n/i) 
*/

void xxs(){
    notpm[1]=1;
    mu[1]=1;
    for(int i=2;i<=maxn;i++){
        if(!notpm[i]){
            pm[++pcnt]=i;mu[i]=-1;
        }
        for(int j=1;j<=pcnt&&pm[j]*i<=maxn;j++){
            notpm[pm[j]*i]=1;
            if(i%pm[j]==0){
                mu[i*pm[j]]=0;
                break;
            }
            mu[i*pm[j]]=-mu[i];
        } 
        mu[i]+=mu[i-1]; 
    }
} 

LL solve(int x,int y){
    x/=k;y/=k;  // 推导公式需要 
    LL tmp=0;
    for(int i=1,j;i<=min(x,y);i=j+1){
        j=min(x/(x/i),y/(y/i)); // 详见上 
        tmp+=1ll*(mu[j]-mu[i-1])*(x/i)*(y/i);
    }
    return tmp;
}

int main(){
//  freopen("2301.in","r",stdin);
//  freopen("2301.out","w",stdout);
    scanf("%d",&tes);
    xxs(); 
    while(tes--){
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        printf("%lld\n",solve(b,d)-solve(b,c-1)-solve(a-1,d)+solve(a-1,c-1));
    }

    return 0;
}
发布了88 篇原创文章 · 获赞 16 · 访问量 1万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 编程工作室 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览