【bzoj3529】[Sdoi2014]数表 莫比乌斯反演+离线+树状数组

题目描述

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

输入

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

输出

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

样例输入

2
4 4 3
10 10 5

样例输出

20
148


题解

莫比乌斯反演+离线+树状数组

设 $f(n)$ 表示 $n$ 的约数和,当约数和大于 $a$ 时将其视作0。

那么我们要求的就是 $sumlimits_{i=1}^nsumlimits_{j=1}^mf(gcd(i,j))$ 

于是有以下公式推导(为了方便,以下默认 $nle m$):

$ sumlimits_{i=1}^nsumlimits_{j=1}^mf(gcd(i,j))\=sumlimits_{d=1}^nf(d)sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=d]\=sumlimits_{d=1}^nf(d)sumlimits_{i=1}^{lfloorfrac nd floor}sumlimits_{j=1}^{lfloorfrac md floor}[gcd(i,j)=1]\=sumlimits_{d=1}^nf(d)sumlimits_{i=1}^{lfloorfrac nd floor}sumlimits_{j=1}^{lfloorfrac md floor}sumlimits_{p|i&p|j}mu(p)\=sumlimits_{d=1}^nf(d)sumlimits_{p=1}^{lfloorfrac nd floor}mu(p)lfloorfrac n{dp} floorlfloorfrac m{dp} floor$

然后再设 $D=dp$ 则有:

$ sumlimits_{d=1}^nf(d)sumlimits_{p=1}^{lfloorfrac nd floor}mu(p)lfloorfrac n{dp} floorlfloorfrac m{dp} floor\=sumlimits_{D=1}^nlfloorfrac nD floorlfloorfrac mD floorsumlimits_{d|D}f(d)mu(frac Dd)$

此时如果所有的询问的 $a$ 都相同,这道题就做完了。预处理出 $f(n)$ 和 $mu(n)$ ,并用 $O(nln n)$ 的时间内算出所有的 $sumlimits_{d|n}f(d)mu(frac nd)$ ,求出其前缀和。然后对每个询问枚举商,分块处理即可。

那么如果 $a$ 不同呢?考虑离线,把所有询问按照 $a$ 从小到大排序,把所有数按照 $f$ 值从小到大排序,对于每一次 $a$ 的改变,在所有数中移动指针 $p$。那么每次指针 $p$ 移动就相当于是把 $f(p)$ 的贡献都算上。由于每个数只被扫到一次,因此只需要修改 $nln n$ 个 $sumlimits_{d|n}f(d)mu(frac nd)$ 。

而我们需要动态维护前缀和,需要支持修改,因此使用树状数组来维护即可。

模 $2^{30}$ 的话可以先unsigned int自然溢出(相当于模$2^{31}$),最后答案再对 $2^{30}-1$取与。

时间复杂度 $O(nlog^2n+Qsqrt nlog n)$

#include <cstdio>
#include <algorithm>
#define N 100010
#define n 100000
using namespace std;
typedef unsigned int ui;
struct data
{
	ui p , q , v , pos;
	bool operator<(const data &a)const {return v < a.v;}
}a[N];
ui mu[N] , prime[N] , tot , np[N] , s[N] , id[N] , f[N] , ans[N];
bool cmp(ui a , ui b)
{
	return s[a] < s[b];
}
inline void add(ui x , ui a)
{
	ui i;
	for(i = x ; i <= n ; i += i & -i) f[i] += a;
}
inline ui query(ui x)
{
	ui i , ans = 0;
	for(i = x ; i ; i -= i & -i) ans += f[i];
	return ans;
}
int main()
{
	ui m , i , j , last , p = 1;
	int x;
	mu[1] = 1;
	for(i = 2 ; i <= n ; i ++ )
	{
		if(!np[i]) mu[i] = -1 , prime[++tot] = i;
		for(j = 1 ; j <= tot && i * prime[j] <= n ; j ++ )
		{
			np[i * prime[j]] = 1;
			if(!(i % prime[j]))
			{
				mu[i * prime[j]] = 0;
				break;
			}
			else mu[i * prime[j]] = -mu[i];
		}
	}
	for(i = 1 ; i <= n ; i ++ )
	{
		id[i] = i;
		for(j = i ; j <= n ; j += i)
			s[j] += i;
	}
	sort(id + 1 , id + n + 1 , cmp);
	scanf("%u" , &m);
	for(i = 1 ; i <= m ; i ++ )
	{
		scanf("%u%u%u" , &a[i].p , &a[i].q , &x) , a[i].pos = i;
		if(x > 0) a[i].v = x;
		else a[i].v = 0;
	}
	sort(a + 1 , a + m + 1);
	for(i = 1 ; i <= m ; i ++ )
	{
		while(p <= n && s[id[p]] <= a[i].v)
		{
			for(j = 1 ; j * id[p] <= n ; j ++ )
				add(j * id[p] , mu[j] * s[id[p]]);
			p ++ ;
		}
		for(j = 1 ; j <= a[i].p && j <= a[i].q ; j = last + 1)
			last = min(a[i].p / (a[i].p / j) , a[i].q / (a[i].q / j)) , ans[a[i].pos] += (a[i].p / j) * (a[i].q / j) * (query(last) - query(j - 1));
	}
	for(i = 1 ; i <= m ; i ++ ) printf("%u
" , ans[i] & 0x7fffffff);
	return 0;
}

 

原文地址:https://www.cnblogs.com/GXZlegend/p/7889121.html