1. 程式人生 > 其它 >演算法提高課

演算法提高課

達達正在破解一段密碼,他需要回答很多類似的問題:

對於給定的整數 a,ba,b 和 dd,有多少正整數對 x,yx,y,滿足 x≤a,y≤bx≤a,y≤b,並且 gcd(x,y)=dgcd(x,y)=d。

作為達達的同學,達達希望得到你的幫助。

輸入格式

第一行包含一個正整數 nn,表示一共有 nn 組詢問。

接下來 nn 行,每行表示一個詢問,每行三個正整數,分別為 a,b,da,b,d。

輸出格式

對於每組詢問,輸出一個正整數,表示滿足條件的整數對數。

資料範圍

1≤n≤500001≤n≤50000,
1≤d≤a,b≤500001≤d≤a,b≤50000

輸入樣例:

2
4 5 2
6 4 3

輸出樣例:

3
2

提示:gcd(x,y)gcd(x,y) 返回 x,yx,y 的最大公約數。

核心思路:首先我們需要知道莫比烏斯函式的相關定義:

我們對於原式子化簡可得:x'<=\(\frac{a}{d}\),y‘<=\(\frac{b}{d}\),gcd(x',y')=1;所以整個問題轉換為了求x'和y'互質的對數,在一定的範圍內。這裡可以用到容斥原理了,也就是合法的對數=總對數-不合法的對數。

總對數很好求:我們令a'=\(\frac{a}{d}\),b'=\(\frac{b}{d}\),所以總對數就是a'b';

那不合法的對數我們怎麼轉換為呢,我們就分情況:找出一個公共的質因子,找出兩個公共的質因子......

使用數學符號就可以表示為:\(\frac{a'}{2}*\frac{b'}{2}\),\(\frac{a'}{6}\frac{b'}{6}\)......。然後我們需要注意下,我們使用容斥原理是需要找集合,所以我們可以吧把只有一個質因子的定義為\(S_1\),有兩個定義為\(S_2\),以此類推。

所以我們使用容斥原理的總表示就是:

其實我們會發現這個是正好符合我們莫比烏斯函式的定義的。因為我們規定只要有一個質因子就是-1(次數不可以大於1),兩個就是1,以此類推。就是\((-1)^k\).
接下來就是數論分塊的思想了。不記得一定要去看部落格,順便看下莫比烏斯反演的部落格,加深印象.
因為我們的\(\frac{a'}{i}\)

\(\frac{b'}{i}\)在一定的區間內這個函式值是不會變的,所以我們需要求一下字首和u(x)函式的。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 6e5 + 10;
int primes[N], st[N];
int mobius[N],sum[N];
int cnt;
void Init()
{
    mobius[1]=1;
	for (int i = 2;i < N;i++)
	{
		if (!st[i])
		{
			primes[cnt++] = i;
			mobius[i] = -1;
		}
		for (int j = 0;primes[j] * i < N;j++)
		{
			int t = primes[j] * i;
			st[t] = 1;
			if (i % primes[j] == 0)//說明至少含有兩個primes[j]
			{
				mobius[t] = 0;
				break;
			}
			mobius[t] = -mobius[i];
		}
	}
	for (int i = 1;i < N;i++)
		sum[i] = sum[i - 1] + mobius[i];
}
int main()
{
	Init();
	int T;
	cin >> T;
	while (T--)
	{
		LL res = 0;
		int a, b, d;
		cin >> a >> b >> d;
		a /= d;
		b /= d;
		int n = min(a, b);//這個和我們的定義有關可以看下那個推導的公式
		for (int l = 1, r;l <= n;l = r + 1)
		{
			r = min(n, min(a / (a / l), b / (b / l)));//這裡有人差不多,反正可以理解為取交集
			res += (sum[r] - sum[l - 1]) * (LL)(a / l) * (b / l);
		}
		cout << res << endl;
	}
}