bzoj2301 莫比乌斯函数

http://www.lydsy.com/JudgeOnline/problem.php?id=2301

其实就是求：

b     d

∑     ∑    [gcd(x,y)==k]

x=a y=c

然后：

b   d
∑  ∑  gcd(x,y)==k
x=a y=c
            a   b
令f[a][b]=  ∑  ∑  gcd(x,y)=k
            x=1 y=1 


ans=f[b][d]+f[a-1][c-1]-f[b][c-1]-f[a-1][d]
             a   b
问题转化为求 ∑  ∑  [gcd(x,y)==d ]
             x=1 y=1
  a'  b'
求∑  ∑   [gcd(x,y)==1]        a'=a/d,b'=b/d 
  x=1 y=1
 
 a'  b'
 ∑  ∑    ∑   u(d)
x=1  y=1  d|gcd(x,y)




min(a',b')
 ∑      u(d)* [a'/d]*[b'/d]

d=1 

做法：预处理线性筛除莫比乌斯函数，然后采用分块。

可以证明[a'/d]的个数是有限的，同理b'/d的个数也是有限的，并且他们在一定范围内是连续相等的，对于连续相等的块，我们可以直接做出来。

对于分块，我们假设当前是从满足【a'/d】=i的第一个d开始的，找到最后一个【a'/d】=i  可以 a'/(a'/d)  分母是整数部分，相除得到d的最大值，从d到这个最大值这一块都是连续相等的a'/d

同理，我们也要找到连续相等的b'/d，把这两个都连续相等的一段直接做出来，由于a'/d不会超过2sqrt(a')个数，复杂度就大大优化了。

（证明：首先[1,sqrt(a)]中,a'/d的值不超过sqrt(a)个；  同理，在[sqrt(a),a]中， 每个因子都会唯一对应一个在sqrt(a)中的因子，所以他也不会超过sqrt(a)个，加起来不超过2*sqrt(a)个 ）

本蒟蒻不明白为什么有大神跑得飞快，我的跑了10s多~~~~（排160）

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<queue>
#include<cmath>
using namespace std;
const int maxn=50000+20;
/*
b   d
∑  ∑  gcd(x,y)==k
x=a y=c
            a   b
令f[a][b]=  ∑  ∑  gcd(x,y)=k
            x=1 y=1 

ans=f[b][d]+f[a-1][c-1]-f[b][c-1]-f[a-1][d]
             a   b
问题转化为求 ∑  ∑  [gcd(x,y)==d ]
             x=1 y=1
  a'  b'
求∑  ∑   [gcd(x,y)==1]        a'=a/d,b'=b/d 
  x=1 y=1
 
 a'  b'
 ∑  ∑    ∑   u(d)
x=1  y=1  d|gcd(x,y)


min(a',b')
 ∑      u(d)* [a'/d]*[b'/d]
d=1 
*/
int u[maxn];
bool flag[maxn];
int prime[maxn];
int sum[maxn];
int tot;
void init()
{
	tot=0;
	memset(flag,0,sizeof(flag));
	memset(sum,0,sizeof(sum));
	u[1]=1;
	for(int i=2;i<=50000;i++)
	{
		if(!flag[i])
		{
			prime[++tot]=i;
			u[i]=-1;
		}
		for(int j=1;j<=tot&&i*prime[j]<=50000;j++)
		{
			flag[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				u[i*prime[j]]=0;
				break;
			}
			else u[i*prime[j]]=u[i]*u[prime[j]];
		}
	}
	for(int i=1;i<=50000;i++)sum[i]=sum[i-1]+u[i];
}
int k;
int work(int a,int b)
{
	a/=k,b/=k;
	int m=min(a,b);
	int pos;
	int ans=0;
	for(int i=1;i<=m;i=pos+1)
	{
		pos=min(a/(a/i),b/(b/i));
		ans+=(sum[pos]-sum[i-1])*(a/i)*(b/i);
	}
	return ans;
}
int n,a,b,c,d;
int main()
{
	init();
	scanf("%d",&n);
	while(n--)
	{
		scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
		int ans=work(b,d)+work(a-1,c-1)-work(a-1,d)-work(b,c-1);
		printf("%d\n",ans);
	}
	return 0;
}