1. 程式人生 > 其它 >【ybtoj高效進階 21177】小小網格(杜教篩)(數論分塊)(莫比烏斯反演)

【ybtoj高效進階 21177】小小網格(杜教篩)(數論分塊)(莫比烏斯反演)

給你求 ∑i=1\~n∑j=1\~mφ(gcd(i,j))。

小小網格

題目連結:ybtoj高效進階 21177

題目大意

給你求 ∑i=1~n∑j=1~mφ(gcd(i,j))。

思路

重新寫好看點:
\(\sum\limits_{i=1}^n\sum\limits_{j=1}^m\varphi(\gcd(i,j))\)
\(\sum\limits_{d}\varphi(d)\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=d]\)
\(\sum\limits_{d}\varphi(d)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]\)


\(\sum\limits_{d}\varphi(d)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum\limits_{p|gcd(i,j)}\mu(p)\)
\(\sum\limits_{d}\varphi(d)\sum\limits_{p}\sum\limits_{i=1,i|p}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1,j|p}^{\left\lfloor\frac{m}{d}\right\rfloor}\mu(p)\)

\(\sum\limits_{d}\varphi(d)\sum\limits_{p}\mu(p)\cdot\left\lfloor\dfrac{n}{dp}\right\rfloor\cdot\left\lfloor\dfrac{m}{dp}\right\rfloor\)

然後設 \(T=pd\),然後有:
\(\sum\limits_{d}\varphi(d)\sum\limits_{p}\mu(p)\cdot\left\lfloor\dfrac{n}{T}\right\rfloor\cdot\left\lfloor\dfrac{m}{T}\right\rfloor\)
\(\sum\limits_{T=1}^{\min(n,m)}\left\lfloor\dfrac{n}{T}\right\rfloor\cdot\left\lfloor\dfrac{m}{T}\right\rfloor\cdot(\sum\limits_{d|T}\varphi(d)\mu(\dfrac{T}{d}))\)

然後你會發現右邊就是一個狄利克雷卷積。
那我們設 \(g=\varphi*\mu\)(積性函式,因為是兩個積性函式相乘),然後式子就是:
\(\sum\limits_{T=1}^{\min(n,m)}\left\lfloor\dfrac{n}{T}\right\rfloor\cdot\left\lfloor\dfrac{m}{T}\right\rfloor\cdot g(d)\)

然後你發現左邊兩個可以數論分塊,然後你考慮右邊的能不能求字首和。
看到是積性函式,考慮用杜教篩,設 \(S(n)=\sum\limits_{i=1}^ng(i)\)

然後你就要找一個合適的 \(f\),不難看出可以用 \(I*I=d\)
然後 \(g*d=\varphi*\mu*I*I=(\varphi*I)*(\mu*I)=id*\epsilon=id\)

然後就有 \(d(1)S(n)=\sum\limits_{i=1}^nid(i)-\sum\limits_{i=2}^nd(i)S(\left\lfloor\dfrac{n}{i}\right\rfloor)=\dfrac{n(n+1)}{2}-\sum\limits_{i=2}^nd(i)S(\left\lfloor\dfrac{n}{i}\right\rfloor)\)

然後 \(d(i)\) 的字首和可以數論分塊,小於 \(n^{\frac{2}{3}}\) 的預處理字首和,大於的直接根號求。

然後小小卡個常即可。

程式碼

#include<map>
#include<cstdio>
#define ll long long
#define mo 1000000007

using namespace std;

const int Maxn = 2000001;
int n, m;
ll d[Maxn + 1], g[Maxn + 1], ans;
int low[Maxn + 1], phi[Maxn + 1], prime[Maxn + 1];
map <int, ll> ans_g;

void init() {//杜教篩的預處理
	g[1] = 1; phi[1] = 1; low[1] = 1;
	for (int i = 2; i <= Maxn; i++) {
		if (!low[i]) phi[i] = i - 1, low[i] = i, prime[++prime[0]] = i, g[i] = i - 2;
		for (int j = 1; j <= prime[0] && i * prime[j] <= Maxn; j++) {
			if (i % prime[j]) low[i * prime[j]] = prime[j], phi[i * prime[j]] = phi[i] * (prime[j] - 1), g[i * prime[j]] = g[i] * g[prime[j]];
				else {
					low[i * prime[j]] = low[i] * prime[j], phi[i * prime[j]] = phi[i] * prime[j];
					if (low[i] != i) g[i * prime[j]] = g[i / low[i]] * g[low[i] * prime[j]];
						else g[i * prime[j]] = phi[i * prime[j]] + phi[i] * (-1);
					break;
				}
		}
	}
	for (int i = 1; i <= Maxn; i++)
		for (int j = 1; i * j <= Maxn; j++)
			d[i * j]++;
	for (int i = 1; i <= Maxn; i++)
		d[i] += d[i - 1], g[i] += g[i - 1];
}

ll sum_d(ll n) {//數論分塊
	if (n <= Maxn) return d[n];
	ll re = 0;
	for (ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		re += (r - l + 1) * (n / l);
	}
	return re;
}

ll sum_g(ll n) {//杜教篩
	if (n <= Maxn) return g[n];
	if (ans_g[n]) return ans_g[n];
	ll re = n * (n + 1) / 2;
	for (ll l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		re -= (sum_d(r) - sum_d(l - 1)) * sum_g(n / l);
	}
	return ans_g[n] = re;
}

int main() {
//	freopen("mesh.in", "r", stdin);
//	freopen("mesh.out", "w", stdout);
	
	scanf("%d %d", &n, &m);
	
	init();
	
	for (int l = 1, r; l <= n && l <= m; l = r + 1) {//數論分塊
		r = min(n / (n / l), m / (m / l));
		ans = (ans + 1ll * (n / l) * (m / l) * (sum_g(r) - sum_g(l - 1))) % mo;
	}
	
	printf("%lld", ans);
	
	return 0;
}