URAL - 1907 Coffee and Buns 容斥原理+数论

http://acm.hust.edu.cn/vjudge/problem/viewProblem.action?id=33020

将题目抽象为

输入 a n 判断 从 b = 1 ~ n 中 满足 gcd( ( a * a + b * b ) , 4 * a + 4 * b ) > 1 的个数。

下面进行推导：

若 a b 奇偶性相同，显然满足。gcd是2的倍数。 

若a b  奇偶性不同，

开始推导：

因为 a * a + b * b 是奇数，

所以

gcd( a * a + b * b , 4 * a + 4 * b ) > 1 等价于

gcd( a * a + b * b , a + b ) > 1 等价于

根据整除性质，若存在 e > 1 ，满足 e | ( a + b ) ， 那么 e | ( a + b ) ^ 2 反之亦然

则等价于

gcd( a * a + b * b , a * a + b * b + 2 * a * b ) > 1 等价于

gcd( 2 * a * a + 2 * b * b , a * a + b * b + 2 * a * b ) > 1等价于

gcd( a * a + b * b - 2 * a * b , a * a + b * b + 2 * a * b ) > 1 等价于

gcd( a - b , a + b ) > 1 等价于

gcd( 2 * a , a + b ) > 1

因为 a+b是奇数

等价于

gcd( a , a +b ) > 1等价于

gcd( a , b ) > 1

推导结束

下面对a分解质因数，找出n个数中不同奇偶性且满足gcd( a , b ) > 1的个数

最后加上 相同奇偶性的 个数 就是答案。

#include <iostream>
#include <cmath>
#define rise(i,a,b,c) for(int i=a;i<=b;i+=c)
long long gcd( long long a , long long b )
{
	if( a % b == 0 )
    return b;
    return gcd( b , a % b );
}
using namespace std;
int main()
{
	long long a, n;
	long long p[200] , pk = 0;
	cin >> a >> n;
	if( n == 1 )
	{
		cout << a % 2;
		return 0;
	}
	if( a == 1 )
	{
		cout << ( n + 1 ) / 2;
		return 0;
	}
	bool qi = true;
	if( a % 2 == 0 )
	{
		qi = false;
		while( a % 2 == 0 )
		a /= 2;
	}
	rise( i , 3 , sqrt(a) , 2 )
	{
		if( a == 1 ) break;
		if( a % i == 0 )
		{
			p[++pk] = i;
			if( qi ) p[pk] *= 2;
			while( a % i == 0 )
			a /= i;
		}
	}
	if( a > 1 ) 
	{
		p[++pk] = a;
		if( qi ) p[pk] *= 2;
	}

	long long ans =0;
	rise( cnt , 1 , (1<<pk)-1 , 1 )
	{
		long long j = 1;
		long long cnt1 = 0;
		rise( i , 0 , pk - 1 , 1 )
		{
			if( cnt&(1<<i) )
			{
				
				j = j / gcd( j , p[i+1] ) * p[i+1];
				cnt1 ++;
				if( j > n ) break;
			}
		}
		if( j > n ) continue;
		if( qi == 0 )
		{
			if( cnt1 & 1 )
			ans += ( n / j + 1 ) / 2;
			else ans = ans - ( n / j + 1 ) / 2;  	
		} 
		else
		{
			if( cnt1 & 1 )
			ans += ( n / j );
			else ans = ans - ( n / j );  
		}
	}
	if( qi == 0 ) cout << ans + n / 2 << endl;
	else cout << ans + ( n - n / 2 ) << endl;
	return 0;
}