首页 > 其他 > 详细

HDU1695 GCD 数论之 莫比乌斯反演

时间:2014-03-01 15:53:43      阅读:537      评论:0      收藏:0      [点我收藏+]

做了一段时间的DP,继续回头啃数论了,这是一道莫比乌斯反演的题目,比较繁琐

给你a,b,c,d,k五个数,其中a=c=1固定的,让你从[a,b]中找出x,[c,d]中找出y,是的gcd(x,y) == k,注意gcd(x,y) 与gcd(y,x)归为同一种,问一共能找到多少组x,y;


分析:


因为gcd(x,y) = k,充要条件gcd(x/k,y/k) = 1,所以区间可以缩小为[1/k,b/k],[1/k,d/k],注意此处的d是题目中的d,不是数论中默认的d最大为公约数含义,这样酒吧问题转换为 去两个区间中的元素x,y使得gcd(x,y) = 1,因为gcd(y,x)和gcd(x,y)是相同的。假设x<y,可以先取小区间进行操作,由于gcd(x,y) = 1,那么可以利用欧拉函数 区间计算为 ans+=phi[1]+phi[2]+……+phi[b]

对于[b/k+1,d/k],设y为区间的一个元素,对y进行素因子分解,得到集合{p1,p2,p3……}pi为素数,求的是gcd(x,y)=1的组合,但是反过来可以求cgd(x,y)!=1的组合,这样就是利用了容斥原理进行统计能被这些素数整出的数的个数,最后相减,求补数即可,然后再加上欧拉函数得到的ans值 就是最后的答案了,



#include<iostream>
#include<cstdio>
#include<list>
#include<algorithm>
#include<cstring>
#include<string>
#include<queue>
#include<stack>
#include<map>
#include<vector>
#include<cmath>
#include<memory.h>
#include<set>

#define ll long long

#define eps 1e-8

#define inf 0xfffffff
//const ll INF = 1ll<<61;

using namespace std;

//vector<pair<int,int> > G;
//typedef pair<int,int > P;
//vector<pair<int,int> > ::iterator iter;
//
//map<ll,int >mp;
//map<ll,int >::iterator p;

//#define IN freopen("c:\\Users\\linzuojun\\desktop\\input.txt", "r", stdin)  
//#define OUT freopen("c:\\Users\\linzuojun\\desktop\\output.txt", "w", stdout) 




//筛选法 同时求出了n以内的素数
const int N = 100000 + 5;
int num[N],cnt;
int a,b,c,d,k,sum;

int phi[N],prime[N];
bool vis[N];

int len;

void init()
{
	ll i,j;
	memset(vis,false,sizeof(vis));
	vis[0]=vis[1]=true;
	len = 0;
	for(i=2;i<=N;i++){
		if(!vis[i])
			prime[len++] = i;
		for(int j=0;j<len && i*prime[j] < N;j++) {
			vis[i*prime[j]] = true;
			if(!(i%prime[j]))break;
		}
	} //这段求出了N内的所有素数
	for(i=1;i<=N;i++) {
		phi[i]=i;
	}
	//sum[0]=0;
	//sum[1]=phi[1];
	memset(phi, 0, sizeof(phi));
	phi[1] = 1;
	for(int i = 2; i < N; i++) if(!phi[i]){
		for(int j = i; j < N; j+=i){
			if(!phi[j]) phi[j] = j;
			phi[j] = phi[j] / i * (i-1);
		}
	}
}//这段求出了N内所有数的欧拉函数值


void dfs(int i,int u,int x,int v) {
	if(u == x) {
		sum += b/v;
		return;
	}
	if(i == cnt)return;
	dfs(i+1,u+1,x,v*num[i]);
	dfs(i+1,u,x,v);
	return;
}

int cal() {//容斥原理
	int s = 0;
	for(int i=1;i<=cnt;i++) {
		sum = 0;
		dfs(0,0,i,1);
		if(i&1)
			s += sum;//奇数加
		else
			s -= sum;//偶数减
	}
	return b - s;
}

int main() {
	init();
	int t;
	int Case = 0;
	scanf("%d",&t);
	while(t--) {
		scanf("%d %d %d %d %d",&a,&b,&c,&d,&k);
		if(k == 0) {
			printf("Case %d: %d\n",++Case,0);
			continue;
		}
		b /= k;
		d /= k;
		if(b > d) {
			a = b;
			b = d;
			d = a;
		}
		ll ans = 0;
		for(int i=1;i<=b;i++)
			ans += phi[i];
		for(int i=b+1;i<=d;i++) {
			cnt = 0;
			a = i;
			for(int j=0;j<len && prime[j]*prime[j] <= a;j++)
				if(!(a%prime[j])) {
					num[cnt++] = prime[j];
					do
					a /= prime[j];
					while(a%prime[j] == 0);
				}
				if(a > 1)
					num[cnt++] = a;
				ans += cal();
		}
		printf("Case %d: %I64d\n",++Case,ans);
	}
	return EXIT_SUCCESS;
}


HDU1695 GCD 数论之 莫比乌斯反演,布布扣,bubuko.com

HDU1695 GCD 数论之 莫比乌斯反演

原文:http://blog.csdn.net/yitiaodacaidog/article/details/20158461

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!