1. 程式人生 > >[bzoj3529][莫比烏斯反演][樹狀陣列]數表

[bzoj3529][莫比烏斯反演][樹狀陣列]數表

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 ( g c d

( i , j ) ) \sum F(gcd(i,j)) ,其中 F
( i ) F(i)
表示 i i 的約數和
列舉一下gcd可以知道是 d i j [ g c d ( i , j ) = d ] F ( d ) \sum_d\sum_i\sum_j [gcd(i,j)=d]F(d)
F ( d ) F(d) 提到前面然後後面就是個計數
反演一下可以知道是
d m i n ( n , m ) i m i n ( n , m ) d F ( d ) μ ( i ) n d i m d i \sum_d^{min(n,m)}\sum_{i}^{\frac{min(n,m)}{d}}F(d)*\mu(i)*\lfloor\frac{n}{d*i}\rfloor*\lfloor\frac{m}{d*i}\rfloor
然後就是套路列舉 d i d*i ,計算前面 F ( d ) μ ( i ) \sum F(d)*\mu (i)
加上a的限制的話就離線一下…每次把約數和在a以內的加入 d i d*i
樹狀陣列維護一下就好了…

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<queue>
#include<vector>
#include<ctime>
#include<map>
#include<bitset>
#include<set>
#define LL long long
#define mp(x,y) make_pair(x,y)
#define pll pair<long long,long long>
#define pii pair<int,int>
using namespace std;
inline int read()
{
	int f=1,x=0;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
int stack[20];
inline void write(LL x)
{
	if(x<0){putchar('-');x=-x;}
    if(!x){putchar('0');return;}
    int top=0;
    while(x)stack[++top]=x%10,x/=10;
    while(top)putchar(stack[top--]+'0');
}
inline void pr1(int x){write(x);putchar(' ');}
inline void pr2(LL x){write(x);putchar('\n');}
const int MAXN=100005;
const int N=100000;
LL bit[MAXN];
int lowbit(int x){return x&-x;}
void change(int x,LL c){for(;x<=N;x+=lowbit(x))bit[x]+=c;}
LL qry(int x){LL ret=0;for(;x>=1;x-=lowbit(x))ret+=bit[x];return ret;}

struct md{long long sum;int p;}w[MAXN];
bool cmp(md n1,md n2){return n1.sum<n2.sum;}

int mu[MAXN],pr[MAXN],plen;
bool vis[MAXN];
void getmu()
{
	vis[1]=true;mu[1]=1;
	for(int i=2;i<=N;i++)
	{
		if(!vis[i])mu[i]=-1,pr[++plen]=i;
		for(int j=1;i*pr[j]<=N&&j<=plen;j++)
		{
			vis[i*pr[j]]=true;
			if(!(i%pr[j])){mu[i*pr[j]]=0;break;}
			else mu[i*pr[j]]=-mu[i];
		}
	}
}
struct ask{int n,m,a,opt;}A[MAXN];
bool cmp1(ask n1,ask n2){return n1.a<n2.a;}
LL answer[MAXN];
int main()
{
//	freopen("2.in","r",stdin);
//	freopen("a.out","w",stdout);
	getmu();
	for(int i=1;i<=N;i++)for(int j=1;j*i<=N;j++)w[i*j].sum+=i;
	for(int i=1;i<=N;i++)w[i].p=i;
	sort(w+1,w+1+N,cmp);
	int T=read();for(int i=1;i<=T;i++)A[i].n=read(),A[i].m=read(),A[i].a=read(),A[i].opt=i;
	sort(A+1,A+1+T,cmp1);
	int LA=0;
	for(int i=1;i<=T;i++)
	{
		while(LA<N&&w[LA+1].sum<=A[i].a)
		{
			LA++;
			for(int j=1;j*w[LA].p<=N;j++)change(w[LA].p*j,w[LA].sum*mu[j]);
		}
		int n=A[i].n,m=A[i].m;
		LL ans=0;int lim=min(n,m),nxt;
		for(int j=1;j<=lim;j=nxt+1)
		{
			nxt=min(n/(n/j),m/(m/j));
			ans+=(LL)(n/j)*(m/j)*(qry(nxt)-qry(j-1));
		}
		answer[A[i].opt]=ans%(2147483648LL);
	}
	for(int i=1;i<=T;i++)pr2(answer[i]);
	return 0; 	 	
}