[HDU4947]GCD Array(莫比乌斯反演)

题面

http://acm.hdu.edu.cn/showproblem.php?pid=4947

前置知识

题解

我们可以维护数组f,使得f时刻等于(a{ imes}{mu})

考虑每次操作1 n d v对于数组a的改变量({Delta}a(x)),发现

[{Delta}a(x)=v*[(x,n)=d] ]

[=v[d|x][d|n][({frac{x}{d}},{frac{n}{d}})=1] ]

如果(d{ mid}n),那么这个操作可以直接忽略。所以不妨设d|n,

[{Delta}a(x)=v[d|x][({frac{x}{d}},{frac{n}{d}})=1] ]

利用消gcd的技巧,反演得

[{Delta}a(x)=v[d|x]{sumlimits_{p|{frac{x}{d}},p|{frac{n}{d}}}}{mu}(p) ]

[={sumlimits_{p|{frac{n}{d}}}}v{mu}(p)[d|x][p|{frac{x}{d}}] ]

[={sumlimits_{p|{frac{n}{d}}}}v{mu}(p)[pd|x] ]

此时的({Delta}f)也应该等于({Delta}a{ imes}{mu})

[{Delta}f(x)={sumlimits_{y|x}}{Delta}a({frac{x}{y}}){mu(y)} ]

[={sumlimits_{y|x}}{mu}(y){sumlimits_{p|{frac{n}{d}}}}v{mu}(p)[pd|{frac{x}{y}}] ]

[={sumlimits_{p|{frac{n}{d}}}}v{mu(p)}{sumlimits_{y|x}}{mu(y)}[pd|{frac{x}{y}}] ]

[={sumlimits_{p|{frac{n}{d}}}} v{mu(p)} [pd|x] {sumlimits_{y|{frac{x}{pd}}}} {mu}(y) ]

[={sumlimits_{p|{frac{n}{d}}}} v{mu(p)} [pd|x] [{frac{x}{pd}}=1] ]

[={sumlimits_{p|{frac{n}{d}}}}v{mu(p)}[x=pd] ]

由此可见,对于1 n d v的修改操作,只需要遍历({frac{n}{d}})的所有因子p,然后将(f(pd)+=v{mu(p)})即可。

对于查询操作,同样将结果用f表示:

[{sumlimits_{i=1}^{x}}a(i)={sumlimits_{i=1}^{x}}{sumlimits_{p|i}}f(p) ]

[={sumlimits_{p}}{sumlimits_{i=1}^{x}}[p|i]f(p) ]

[={sumlimits_{p}}{lfloor}{frac{x}{p}}{ floor}f(p) ]

将f的前缀和使用树状数组维护,数论分块即可优化到每次询问(O({sqrt{x}}log x))

但是,注意到这里的区间查询是O(logx)而非普通数论分块的O(1),因此此处可以进一步优化:

设置阈值T,当查询的x小于T时使用暴力求和,否则使用数论分块。

则暴力求和的时间复杂度为O(T),数论分块的时间复杂度为(O({frac{x}{T}}log x))。不难发现置(T={sqrt{xlog x}})时总时间复杂度最低。至此,每次询问的时间复杂度降为(O(sqrt{x log x}))

分析总时间,各预处理均在nlogn级别及以下,(n为数组长度),每次修改是(O({sqrt{n}} log n)),查询是(O(sqrt{n log n}))。故总时间复杂度应为(O(CQ{sqrt{n}} log n)),(C为数据组数),且此处的根号来源于因子个数,其常数极小。

代码

#include<bits/stdc++.h>

using namespace std;

#define ll long long
#define rg register
#define N 50000
#define W 200000

inline ll read(){
	ll s = 0,ww = 1;
	char ch = getchar();
	while(ch < '0' || ch > '9'){if(ch == '-')ww = -1;ch = getchar();}
	while('0' <= ch && ch <= '9'){s = 10 * s + ch - '0';ch = getchar();}
	return s * ww;
}

inline void write(ll x){
	if(x < 0)putchar('-'),x = -x;
	if(x > 9)write(x / 10);
	putchar('0' + x % 10);
}

ll n,Q;
ll a[N+5],f[N+5],lg[W+5];
vector<ll>fac[W+5];

inline void prepro(){ //预处理因子和log 
	for(rg ll i = 1;i <= W;i++)
		for(rg ll j = i;j <= W;j += i)fac[j].push_back(i);
	lg[1] = 0;
	for(rg ll i = 2;i <= W;i++)lg[i] = lg[i>>1] + 1;
}

ll pn;
bool isp[W+5];
ll pri[40000+5],mu[W+5];

inline void Eular(){ //线性筛 
	mu[1] = 1;
	for(rg ll i = 2;i <= W;i++)isp[i] = 1;
	for(rg ll i = 2;i <= W;i++){
		if(isp[i])pri[++pn] = i,mu[i] = -1;
		for(rg ll j = 1;i * pri[j] <= W;j++){
			isp[i*pri[j]] = 0;
			if(i % pri[j])mu[i*pri[j]] = -mu[i];
			else{
				mu[i*pri[j]] = 0;
				break;
			}
		}
	}
}

inline ll lowbit(ll x){
	return x & -x;
}

struct BIT{
	ll b[N+5];
	
	inline void reset(){
		memset(b,0,sizeof(b));
	}
	
	inline void ud(ll x,ll dx){
		for(rg ll i = x;i <= n;i += lowbit(i))b[i] += dx;
	}
	
	inline ll query(ll x){
		ll rt = 0;
		for(rg ll i = x;i;i -= lowbit(i))rt += b[i];
		return rt;
	}
	
	inline ll sum(ll l,ll r){
		return query(r) - query(l - 1);
	}
}B;

int main(){
	prepro(); 
	Eular(); 
	n = read(),Q = read();
	ll t = 0;
	while(n || Q){
		printf("Case #%lld:
",++t);
		B.reset();
		memset(a,0,sizeof(a));
		memset(f,0,sizeof(f));
		while(Q--){
			ll opt = read();
			if(opt == 1){
				ll _n = read(),d = read(),v = read();
				if(_n % d)continue;
				for(rg ll i = 0;i < fac[_n/d].size();i++){
					ll p = fac[_n/d][i];
					if(p * d > n)break;
					f[p*d] += mu[p] * v;
					B.ud(p * d,mu[p] * v); 
				}
			}
			else{
				ll x = read();
				ll T = (ll)round(sqrt(x*lg[x])); 
				if(T > x)T = x;
				ll ans = 0;
				for(rg ll i = 1;i <= T;i++)ans += f[i] * (x / i); //sqrt(xlogx)以下暴力求和 
				ll L,R = T;
				while(R != x){ //sqrt(xlogx)以上数论分块 
					L = R + 1;
					R = x / (x / L);
					ans += B.sum(L,R) * (x / L);
				}
				write(ans),putchar('
');
			}
		}
		n = read(),Q = read();	
	}
	return 0;
}
原文地址:https://www.cnblogs.com/xh092113/p/12272719.html