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;
}