BZOJ_3529_[Sdoi2014]数表_莫比乌斯反演+树状数组

Description

有一张 n×m 的数表,其第 i 行第 j 列(1 <= i <= n, 1 <= j <= m)的数值为
能同时整除 i 和 j 的所有自然数之和。给定 a , 计算数表中不大于 a 的数之和。

Input

输入包含多组数据。
输入的第一行一个整数Q表示测试点内的数据组数
接下来Q行,每行三个整数n,m,a(|a| < =10^9)描述一组数据。
1 < =N.m < =10^5  , 1 < =Q < =2×10^4

Output

对每组数据,输出一行一个整数,表示答案模2^31的值。

Sample Input

2
4 4 3
10 10 5

Sample Output

20
148

考虑如果a相等,多次询问怎么做。
设$f(x)=\{^{0,x>a}_{\sigma(x),x\le a}$

$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}f(gcd(i,j))$

$\sum\limits_{d=1}^{n}f(d)\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[gcd(i,j)=d]$

$\sum\limits_{d=1}^{n}f(d)\sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{m/d}[gcd(i,j)=1]$

$\sum\limits_{d=1}^{n}f(d)\sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{m/d}\sum\limits_{p|gcd(i,j)}\mu(p)$


$\sum\limits_{d=1}^{n}f(d)\sum\limits_{p|n}\mu(p)n/dp*m/dp$

$\sum\limits_{Q=1}^{n}n/Q*m/Q\sum\limits_{d|Q}f(d)\mu(Q/d)$

$g(n)=\sum\limits_{d|n}f(d)\mu(n/d)$
 
预处理出g的前缀和就可以在$O(\sqrt n)$的时间内求出。
现在每次a的值不一样,可以把询问按a排个序。
每插入一个a值相当于修改一些g上的值,且总共需要修改$O(nln)$次。
然后每次还要快速的求g的前缀和,树状数组维护一下g的值即可。
 
代码:
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cstdlib>
using namespace std;
#define N 100050
#define M 100000
typedef unsigned int un;
int prime[N],cnt,miu[N],ys[N],ysh[N],id[N];
un c[N],ans[N];
bool vis[N];
struct A {
	int id,n,m,a;
	bool operator < (const A &x) const {
		return a<x.a;
	}
}q[N];
void fix(int x,un v) {for(;x<=M;x+=x&(-x)) c[x]+=v;}
un inq(int x) {un re=0; for(;x;x-=x&(-x)) re+=c[x]; return re;}
inline bool cmp(const int &x,const int &y) {return ysh[x]<ysh[y];}
void init() {
	int i,j; miu[1]=1;
	for(id[1]=1,i=2;i<=M;i++) {
		if(!vis[i]) {
			prime[++cnt]=i; miu[i]=-1;
		}
		for(j=1;j<=cnt&&i*prime[j]<=M;j++) {
			int y=i*prime[j]; vis[y]=1;
			if(i%prime[j]==0) {miu[y]=0;break;}
			miu[y]=-miu[i];
		}
		id[i]=i;
	}
	for(i=1;i<=M;i++) for(j=i;j<=M;j+=i) ysh[j]+=i;
	sort(id+1,id+M+1,cmp);
}
void add(int x) {
	int i;
	for(i=x;i<=M;i+=x) {
		fix(i,ysh[x]*miu[i/x]);
	}
}
un solve(un n,un m) {
	int i,lst; if(n>m) swap(n,m); un re=0;
	for(i=1;i<=n;i=lst+1) {
		lst=min(n/(n/i),m/(m/i));
		re+=(n/i)*(m/i)*(inq(lst)-inq(i-1));
	}
	return re;
}
int main() {
	init();
	int Q;
	scanf("%d",&Q);
	int i,j=1;
	for(i=1;i<=Q;i++) scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a),q[i].id=i;
	sort(q+1,q+Q+1);
	for(i=1;i<=Q;i++) {
		while(ysh[id[j]]<=q[i].a&&j<=M) add(id[j]),j++;
		ans[q[i].id]=solve(q[i].n,q[i].m);
	}
	for(i=1;i<=Q;i++) printf("%u\n",ans[i]&((1u<<31)-1));
}

 

转载于:https://www.cnblogs.com/suika/p/9429704.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值