主頁 > 後端開發 > 高中生都能看懂的莫比烏斯反演

高中生都能看懂的莫比烏斯反演

2020-09-14 09:43:56 後端開發

寫在前面

額,標題好像少了一個不字……應該是高中生都不能看懂的莫比烏斯反演

這是蒟蒻第一次寫這么長的博文

\(gyh\ nb\)\(\text{OI-Wiki}\ nb\)

如果覺得寫得湊合就點個支持吧\(qwq\)

前置知識

積性函式、狄利克雷卷積、數論分塊(這一篇去找gyh吧我講也講不好)(有空慢慢補)

Mobius函式

定義

莫比烏斯函式\(\mu(n)\)定義為:

\[\mu(n)=\left\{ \begin{aligned} 1, & n=1 \\ (-1)^{s}, & n=p_1p_2…p_s(s為n的本質不同的質因子個數)\\0,&n有平方因子\end{aligned} \right.\]

其中\(p_1p_2…,p_s\)是不同素數,

可以看出,當\(n\)沒有平方因子時,\(\mu(n)\)非零,

\(\mu\)也是積性函式,

性質

莫比烏斯函式具有如下的性質:

\[\sum_{d|n}\mu(d)=\epsilon(n)=[n=1] \]

使用狄利克雷卷積來表示,即

\[\mu*1=\epsilon \]

證明:

\(n=1\)時顯然成立,

\(n>1\),設\(n\)\(s\)個不同的素因子,由于\(\mu(d)\neq0\)當且僅當\(d\)無平方因子,故\(d\)中每個素因子的指數只能為\(0\)\(1\)

\(n\)的某個因子\(d\),有\(\mu(d)=(-1)^i\),則它由\(i\)個本質不同的質因子構成,因為質因子的總數為\(s\),所以滿足上式的因子數有\(C_s^i\)個,

所以我們就可以對于原式,轉化為列舉\(\mu(d)\)的值,同時運用二項式定理,故有

\[\sum_{d|n}\mu(d)=\sum_{i=0}^{s}C_s^i\times(-1)^i=\sum_{i=0}^{s}C_s^i\times(-1)^i\times 1^{s-i}=(1+(-1))^s \]

該式在\(s=0\)\(n=1\)時為1,否則為\(0\)

求莫比烏斯函式

因為莫比烏斯函式是積性函式,所以可以用線性篩

int n, cnt, p[A], mu[A];
bool vis[A];

void getmu() {
    mu[1] = 1;
    for (int i = 2; i <= n; i++) {
	if (!vis[i]) mu[i] = -1, p[++cnt] = i;
	for (int j = 1; j <= cnt && i * p[j] <= n; j++) {
	    vis[i * p[j]] = 1;
	    if (i % p[j] == 0) { mu[i * p[j]] = 0; break; }
	    mu[i * p[j]] -= mu[i];
	}
    }
}

莫比烏斯反演公式

\(f(n)\)\(g(n)\)為兩個數論函式,

如果有

\[f(n)=\sum\limits_{d|n}g(d) \]

則有

\[g(n)=\sum\limits_{d|n}\mu(d)f(\frac{n}{d}) \]

證明

  • 法一:對原式做數論變換

    1. \(\sum\limits_{d|n}g(d)\)替換\(f(n)\),即

      \[\sum\limits_{d|n}\mu(d)f(\frac{n}{d})=\sum\limits_{d|n}\mu(d)\sum\limits_{k|\frac nd}g(k) \]

    2. 變換求和順序得

      \[\sum\limits_{k|n}g(k)\sum\limits_{d|\frac n k}\mu(\frac nk) \]

    3. 因為\(\sum\limits_{d|n}\mu(d)=[n=1]\),所以只有在\(\frac{n}{k}=1\)\(n=k\)時才會成立,所以上式等價于

      \[\sum\limits_{d|n}[n=k]g(k)=g(n) \]

    得證

  • 法二:利用卷積

    原問題為:已知\(f=g*1\),證明\(g=f*\mu\)

    轉化:\(f*\mu=g*1*\mu=g*\varepsilon=g\)\(1*\mu=\varepsilon\)

    再次得證= =

小性質

\([\gcd(i,j)=1]\Leftrightarrow\sum\limits_{d|\gcd(i,j)}\mu(d)\)

證明

  • 法一:

    \(n=\gcd(i,j)\),那么等式右邊\(=\sum\limits_{d|n}\mu(d)=[n=1]=[\gcd(i,j)=1]=\)等式左邊

  • 法二:

    利用單位函式暴力拆開:\([\gcd(i,j)=1]=\varepsilon(\gcd(i,j))=\sum\limits_{d|\gcd(i,j)}\mu(d)\)

做題思路&&應用

利用狄利克雷卷積可以解決一系列求和問題,常見做法是使用一個狄利克雷卷積替換求和式中的一部分,然后調換求和順序,最終降低時間復雜度,

經常利用的卷積有\(\mu*1=\epsilon\)\(\text{Id}=\varphi*1\)

還是以題為主吧,但是做的題也會單獨寫題解,畢竟要多水幾篇博客的嘛/huaji

洛谷 P2522 [HAOI2011]Problem b

題目鏈接

我的題解

\(n\)組詢問,每次給出\(a,b,c,d,k\),求\(\sum\limits_{x=a}^{b}\sum\limits_{y=c}^{d}[\gcd(x,y)=k]\)

\(f(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j= 1}^{m}[\gcd(i,j)=k]\)

那么根據容斥原理,題目中的式子就轉化成了\(f(b,d)-f(b, c - 1) - f(a - 1,d) + f(a - 1, c - 1)\)

所以我們接下來的問題就轉化為了如何求\(f\)的值

現在來化簡\(f\)的值

  1. 容易得出原式等價于$$\sum\limits_{i = 1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j = 1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i,j) = 1]$$

  2. 因為\(\epsilon(n) =\sum\limits_{d|n}\mu(d)=[n=1]\),由此可將原式化為

    \[\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum\limits_{d|gcd(i,j)}\mu(d) \]

  3. 改變列舉物件并改變列舉順序,先列舉\(d\),得

    \[\sum\limits_{d=1}^{\min(n,m)}\mu(d)\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}[d|i]\sum\limits_{j=1}^{\lfloor\frac{m}{k}\rfloor}[d|j] \]

    也就是說當\(d|i\)\(d|j\)時,\(d|\gcd(i,j)\)

  4. 易得\(1\sim \lfloor\frac{n}{k}\rfloor\)中一共有\(\lfloor\frac{n}{dk}\rfloor\)\(d\)的倍數,同理\(1\sim \lfloor\frac{m}{k}\rfloor\)中一共有\(\lfloor\frac{m}{dk}\rfloor\)\(d\)的倍數,于是原式化為$$\sum\limits_{d=1}^{\min(n,m)}\mu(d)\lfloor\frac{n}{dk}\rfloor\lfloor\frac{m}{dk}\rfloor$$

此時已經可以\(O(n)\)求解,但是過不了,因為有很多值相同的區間,所以可以用數論分塊來做

先預處理\(\mu\),再用數論分塊,復雜度\(O(n+T\sqrt n)\)

我的代碼每次得分玄學,看評測機心情,建議自己寫

/*
Author:loceaner
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

const int A = 1e6 + 11;
const int B = 1e6 + 11;
const int mod = 1e9 + 7;
const int inf = 0x3f3f3f3f;

inline int read() {
	char c = getchar(); int x = 0, f = 1;
	for( ; !isdigit(c); c = getchar()) if(c == '-') f = -1;
	for( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
	return x * f;
}

int n, a, b, c, d, k, cnt, p[A], mu[A], sum[A];
bool vis[A];

void getmu() {
	int MAX = 50010;
	mu[1] = 1;
	for (int i = 2; i <= MAX; i++) {
		if (!vis[i]) mu[i] = -1, p[++cnt] = i;
		for (int j = 1; j <= cnt && i * p[j] <= MAX; j++) {
			vis[i * p[j]] = true;
			if (i % p[j] == 0) break;
			mu[i * p[j]] -= mu[i];
		}
	}
	for (int i = 1; i <= MAX; i++) sum[i] = sum[i - 1] + mu[i];
}

int work(int x, int y) {
	int ans = 0ll;
	int max = min(x, y);
	for (int l = 1, r; l <= max; l = r + 1) {
		r = min(x / (x / l), y / (y / l));
		ans += (1ll * x / (l * k)) * (1ll * y / (l * k)) * 1ll * (sum[r] - sum[l - 1]); 
	}
	return ans;
}

void solve() {
	a = read(), b = read(), c = read(), d = read(), k = read();
	cout << work(b, d) - work(a - 1, d) - work(b, c - 1) + work(a - 1, c - 1) << '\n';
}

signed main() {
	getmu();
	int T = read();
	while (T--) solve();
	return 0;
}

洛谷 P3455 [POI2007]ZAP-Queries

題目鏈接

我的題解

\(T\)組詢問,每次詢問求

\[\sum\limits_{x=1}^{a}\sum\limits_{y=1}^{b}[\gcd(x,y)=d] \]

因為我不喜歡用\(x、y、a、b、d\),所以一一對應換成\(i、j、n、m、k\)

直接淦式子

\[\begin{align*}&\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=k]\\=& \sum\limits_{i=1}^{\lfloor \frac nk\rfloor}\sum\limits_{j=1}^{\lfloor\frac mk\rfloor}[\gcd(i,j)=1]\\=&\sum\limits_{i=1}^{\lfloor \frac nk\rfloor}\sum\limits_{j=1}^{\lfloor\frac mk\rfloor}\sum\limits_{d|\gcd(i,j)}\mu(d)\\=&\sum\limits_{i=1}^{\lfloor \frac nk\rfloor}\sum\limits_{j=1}^{\lfloor\frac mk\rfloor}\sum\limits_{d|i}\sum _{d|j}\mu(d)\\=&\sum\limits_{d=1}^{\min(\lfloor\frac nk\rfloor,\lfloor\frac mk\rfloor)}\mu(d)\sum\limits_{i=1}^{\lfloor \frac nk\rfloor}\sum\limits_{j=1}^{\lfloor\frac mk\rfloor}\sum\limits_{d|i}\sum\limits_{d|j}1\\=&\sum\limits_{d=1}^{\min(\lfloor\frac nk\rfloor,\lfloor\frac mk\rfloor)}\mu(d)\sum\limits_{i=1}^{\lfloor \frac nk\rfloor}\sum\limits_{d|i}1\sum\limits_{j=1}^{\lfloor\frac mk\rfloor}\sum\limits_{d|j}1\\=&\sum\limits_{d=1}^{\min(\lfloor\frac nk\rfloor,\lfloor\frac mk\rfloor)}\mu(d)\lfloor\frac{\lfloor \frac nk\rfloor}d\rfloor\lfloor\lfloor\frac{\frac mk\rfloor}d\rfloor\\=&\sum\limits_{d=1}^{\min(\lfloor\frac nk\rfloor,\lfloor\frac mk\rfloor)}\mu(d)\lfloor\frac n{kd}\rfloor\lfloor\frac m{kd}\rfloor\end{align*} \]

現在就可以每次詢問\(O(n)\)做這道題了

但是跑不過啊,不過顯然可以數論分塊,所以我們就可以\(O(\sqrt n)\)回答每次詢問了

/*
Author:loceaner
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

const int A = 5e5 + 11;
const int B = 1e6 + 11;
const int mod = 1e9 + 7;
const int inf = 0x3f3f3f3f;

inline int read() {
	char c = getchar(); int x = 0, f = 1;
	for( ; !isdigit(c); c = getchar()) if(c == '-') f = -1;
	for( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
	return x * f;
}

int n, m, k, mu[A], p[A], sum[A], cnt;
bool vis[A];

void getmu(int n) {
	mu[1] = 1;
	for (int i = 2; i <= n; i++) {
		if (!vis[i]) p[++cnt] = i, mu[i] = -1;
		for (int j = 1; j <= cnt && i * p[j] <= n; j++) {
			vis[i * p[j]] = 1;
			if (i % p[j] == 0) break;
			mu[i * p[j]] -= mu[i];
		}
	}
	for (int i = 1; i <= n; i++) sum[i] = sum[i - 1] + mu[i];
}

int solve(int n, int m, int k) {
	int ans = 0, maxn = min(n, m);
	for (int l = 1, r; l <= maxn; l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans += (sum[r] - sum[l - 1]) * (n / (k * l)) * (m / (k * l));
	}
	return ans;
}

int main() {
	getmu(50000);
	int T = read();
	while (T--) {
		n = read(), m = read(), k = read();
		cout << solve(n, m, k) << '\n';
	}
	return 0;
}

洛谷 P1829 [國家集訓隊]Crash的數字表格 / JZPTAB

題目鏈接

我的題解

\[\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\text{lcm}(i,j)(\bmod 20101009) \]

容易想到原式等價于

\[\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{ij}{\gcd(i,j)} \]

列舉 \(i,j\) 的最大公約數 \(d\),顯然 \(\gcd(\frac id,\frac jd)=1\),即 \(\frac id\)\(\frac jd\) 互質

\[\sum\limits_{i=1}^{n}\sum\limits_{j=1}^m\sum\limits_d[d|i][d|j][\gcd(\frac id,\frac jd)=1]\frac{ij}d \]

變換求和順序,改為先列舉 \(d\),則可以得到:

\[\sum\limits_{d=1}\sum\limits_{i=1}^{n}[d|i]\sum\limits_{j=1}^{m}[d|j][\gcd(\frac id,\frac jd)=1]\frac{ij}{d} \]

把后面的 \(d\) 提取出來,對應的兩個 \(\sum\) 改為列舉 \(\frac id,\frac jd\),進一步化簡可以得到:

\[\begin{aligned}&\sum\limits_{d=1}\sum\limits_{i=1}^{n}[d|i]\sum\limits_{j=1}^{m}[d|j][\gcd(\frac id,\frac jd)=1]\frac{ij}{d}\\=&\sum\limits_{d=1}\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]\frac{id\times jd}{d}\\=&\sum\limits_{d=1}\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]ijd\\=&\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]ij\end{aligned} \]

先擱置這個式子,考慮后半部分,

\(sum(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=1]ij\),對其進行化簡,用 \(\varepsilon(\gcd(i,j))\)(即\(\sum\limits_{d|、gcd(i,j)}\mu(d)\)) 替換\([\gcd(i,j)=1]\),得到

\[\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d|\gcd(i,j)}\mu(d)ij \]

轉化為首先列舉約數

\[\sum\limits_{d=1}^{\min(n,m)}\mu(d)\sum\limits_{i=1}^{n}[d|i]\sum\limits_{j=1}^{m}[d|j]ij \]

像上面一樣,后面的式子中考慮把 \(d\) 提取出來,直接列舉滿足條件的數,得到

\[\begin{aligned}&\sum\limits_{d=1}^{\min(n,m)}\mu(d)\sum\limits_{i=1}^{n}[d|i]\sum\limits_{j=1}^{m}[d|j]ij\\&=\sum\limits_{d=1}^{\min(n,m)}\mu(d)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}id\times jd\\&=\sum\limits_{d=1}^{\min(n,m)}\mu(d)\times d^2\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}ij\end{aligned} \]

前半段可以處理前綴和,后半段可以 \(O(1)\) 求,因為設

\[Q(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij=\frac{n\times(n+1)}{2}\times\frac{m\times(m+1)}{2} \]

的話,這個函式顯然可以\(O(1)\)求解

到現在

\[sum(n,m)=\sum\limits_{d=1}^{\min(n,m)}\mu(d)\times d^2\times Q(\lfloor\frac nd \rfloor,\lfloor\frac md\rfloor) \]

可以用數論分塊求解

回帶到原式中

\[\sum\limits_{d=1}^{\min(n, m)}d\times sum(\lfloor\frac nd \rfloor,\lfloor\frac md\rfloor) \]

這樣的話,又可以數論分塊求解了,然后就做完啦

/*
Author:loceaner
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define int long long
using namespace std;

const int A = 1e7 + 11;
const int B = 1e6 + 11;
const int mod = 20101009;
const int inf = 0x3f3f3f3f;

inline int read() {
	char c = getchar(); int x = 0, f = 1;
	for( ; !isdigit(c); c = getchar()) if(c == '-') f = -1;
	for( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
	return x * f;
}

bool vis[A];
int n, m, mu[A], p[B], sum[A], cnt;

void getmu() {
    mu[1] = 1;
    int k = min(n, m);
    for (int i = 2; i <= k; i++) {
        if (!vis[i]) p[++cnt] = i, mu[i] = -1;
        for (int j = 1; j <= cnt && i * p[j] <= k; ++j) {
            vis[i * p[j]] = 1;
            if (i % p[j] == 0) break;
            mu[i * p[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= k; i++) sum[i] = (sum[i - 1] + i * i % mod * mu[i]) % mod;
}

int Sum(int x, int y) { 
	return (x * (x + 1) / 2 % mod) * (y * (y + 1) / 2 % mod) % mod; 
}

int solve2(int x, int y) {
    int res = 0;
    for (int i = 1, j; i <= min(x, y); i = j + 1) {
        j = min(x / (x / i), y / (y / i));
        res = (res + 1LL * (sum[j] - sum[i - 1] + mod) * Sum(x / i, y / i) % mod) % mod;
    }
    return res;
}

int solve(int x, int y) {
    int res = 0;
    for (int i = 1, j; i <= min(x, y); i = j + 1) {
        j = min(x / (x / i), y / (y / i));
        res = (res + 1LL * (j - i + 1) * (i + j) / 2 % mod * solve2(x / i, y / i) % mod) % mod;
    }
    return res;
}

signed main() {
    n = read(), m = read();
    getmu();
    cout << solve(n, m) << '\n';
}

洛谷 P2257 YY的GCD

題目鏈接

給定\(n,m\),求二元組\((x,y)\)的個數,滿足\(1\leq x\leq n,1\leq y\leq m\),且\(gcd(x,y)\)是素數,

\(n,m\leq 10^7\),自帶多組資料,至多\(10^{4}\)組資料,

思路與第一題Problem B類似,在這里不再贅述,只給出代碼= =

#include <cmath> 
#include <cstdio>
#include <cstring> 
#include <iostream>
using namespace std;

const int A = 1e7 + 11; 

inline int read() {
	char c = getchar(); int x = 0, f = 1;
	for ( ; !isdigit(c); c = getchar()) if(c == '-') f = -1; 
	for ( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
	return x * f;
}

bool vis[A];
long long sum[A];
int prim[A], mu[A], g[A], cnt, n, m;

void get_mu(int n) {
    mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) {
            mu[i] = -1;
            prim[++cnt] = i;
        }
        for (int j = 1; j <= cnt && prim[j] * i <= n; j++) {
            vis[i * prim[j]] = 1;
            if (i % prim[j] == 0) break;
            else mu[prim[j] * i] = - mu[i];
        }
    }
    for (int j = 1; j <= cnt; j++)
        for (int i = 1; i * prim[j] <= n; i++) g[i * prim[j]] += mu[i];
    for (int i = 1; i <= n; i++) sum[i] = sum[i - 1] + (long long)g[i];
}

signed main() {
    int t = read();
    get_mu(10000000);
    while (t--) {
        n = read(), m = read();
        if (n > m) swap(n, m);
        long long ans = 0;
        for (int l = 1, r; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            ans += 1ll * (n / l) * (m / l) * (sum[r] - sum[l - 1]);
        }
        cout << ans << '\n';
    }
    return 0;
}

洛谷 P3327 [SDOI2015]約數個數和

題目鏈接

\[\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}d(ij) \]

\(d(x)\)\(x\)的約數個數和

需要用到

\[d(ij)=\sum\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)=1] \]

證明我也不會

然后自己推導吧,在此不再贅述

/*
Author:loceaner
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

const int A = 5e5 + 11;
const int B = 1e6 + 11;
const int mod = 1e9 + 7;
const int inf = 0x3f3f3f3f;

inline int read() {
	char c = getchar(); int x = 0, f = 1;
	for( ; !isdigit(c); c = getchar()) if(c == '-') f = -1;
	for( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
	return x * f;
}

bool vis[A];
int n, m, p[A], mu[A], cnt, sum[A];
long long g[A], ans;

void getmu(int n) {
	mu[1] = 1;
	for (int i = 2; i <= n; i++) {
		if (!vis[i]) mu[i] = -1, p[++cnt] = i;
		for (int j = 1; j <= cnt && i * p[j] <= n; j++) {
			vis[i * p[j]] = 1;
			if (i % p[j] == 0) break;
			mu[i * p[j]] -= mu[i];
		}
	}
	for (int i = 1; i <= n; i++) sum[i] = sum[i - 1] + mu[i];
	for (int i = 1; i <= n; i++) {
		int ans = 0;
		for (int l = 1, r; l <= i; l = r + 1) {
			r = (i / (i / l));
			ans += 1ll * (r - l + 1) * (i / l);
		}
		g[i] = ans;
	}
}

signed main() {
	int T = read();
	getmu(50000);
	while (T--) {
		n = read(), m = read();
		int maxn = min(n, m);
		ans = 0;
		for (int l = 1, r; l <= maxn; l = r + 1) {
			r = min(n / (n / l), m / (m / l));
			ans += 1ll * (sum[r] - sum[l - 1]) * 1ll * g[n / l] * 1ll * g[m / l];
		}
		cout << ans << '\n';
	}
	return 0;
}

洛谷 P4449 于神之怒加強版

\[\sum_{i=1}^{n}\sum_{j=1}^{m}\gcd(i,j)^k(\bmod 1e9+7) \]

還是直接淦式子

\[\begin{align*}&\sum_{i=1}^{n}\sum_{j=1}^{m}\gcd(i,j)^k\\=&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d=1}^{\min(n,m)}d^k[\gcd(i,j)=d]\\=&\sum_{d=1}^{\min(n,m)}d^k\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[\gcd(i,j)=1]\\=&\sum_{d=1}^{\min(n,m)}d^k\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}\sum_{x|\gcd(i,j)}\mu(x)\\=&\sum_{d=1}^{\min(n,m)}d^k\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}\sum_{x|i,x|j}\mu(x)\\=&\sum_{d=1}^{\min(n,m)}d^k\sum_{x=1}^{\min(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor)}\mu(x)\sum_{i=1}^{\lfloor\frac nd\rfloor}[x|i]\sum_{j=1}^{\lfloor\frac md\rfloor}[x|j]\\=&\sum_{d=1}^{\min(n,m)}d^k\sum_{x=1}^{\min(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor)}\mu(x)\lfloor\frac n{dx}\rfloor\lfloor\frac m{dx}\rfloor\end{align*} \]

\(P=dx\),則原式等于

\[\sum_{P=1}^{\min(n,m)}\lfloor\frac n{P}\rfloor\lfloor\frac m{P}\rfloor\sum_{d|P}d^k\mu(\frac Pd) \]

顯然前面的\(\lfloor\frac n{P}\rfloor\lfloor\frac m{P}\rfloor\)部分可以分塊求解,

現在考慮后面的一部分,令

\[g(n)=\sum_{d|n}d^k\mu(\frac nd) \]

容易得出這個函式是積性函式,所以我們就可以線性篩然后求出其前綴和

然后就做完了

/*
Author:loceaner
莫比烏斯反演
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define int long long
using namespace std;

const int A = 5e6 + 11;
const int B = 1e6 + 11;
const int mod = 1e9 + 7;
const int inf = 0x3f3f3f3f;

inline int read() {
	char c = getchar();
	int x = 0, f = 1;
	for( ; !isdigit(c); c = getchar()) if(c == '-') f = -1;
	for( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
	return x * f;
}

bool vis[A];
int T, n, m, k, f[A], g[A], p[A], cnt, sum[A];

inline int power(int a, int b) {
	int res = 1;
	while (b) {
		if (b & 1) res = res * a % mod;
		a = a * a % mod, b >>= 1;
	}
	return res;
}

inline int mo(int x) {
	if(x > mod) x -= mod;
	return x;
}

inline void work() {
	g[1] = 1;
	int maxn = 5e6 + 1;
	for (int i = 2; i <= maxn; i++) {
		if (!vis[i]) { p[++cnt] = i, f[cnt] = power(i, k), g[i] = mo(f[cnt] - 1 + mod); }
		for (int j = 1; j <= cnt && i * p[j] <= maxn; j++) {
			vis[i * p[j]] = 1;
			if (i % p[j] == 0) { g[i * p[j]] = g[i] * 1ll * f[j] % mod; break; }
			g[i * p[j]] = g[i] * 1ll * g[p[j]] % mod;
		}
	}
	for (int i = 2; i <= maxn; i++) g[i] = (g[i - 1] + g[i]) % mod;
}

inline int abss(int x) {
	while (x < 0) x += mod;
	return x;
}

signed main() {
	T = read(), k = read();
	work();
	while (T--) {
		n = read(), m = read();
		int maxn = min(n, m), ans = 0;
		for (int l = 1, r; l <= maxn; l = r + 1) {
			r = min(n / (n / l), m / (m / l));
			(ans += abss(g[r] - g[l - 1]) * 1ll * (n / l) % mod * (m / l) % mod) %= mod;
		}
		ans = (ans % mod + mod) % mod;
		cout << ans << '\n';
	}
	return 0;
}

轉載請註明出處,本文鏈接:https://www.uj5u.com/houduan/33460.html

標籤:C++

上一篇:抽象寵物類的實作 代碼參考

下一篇:STL之deque

標籤雲
其他(157675) Python(38076) JavaScript(25376) Java(17977) C(15215) 區塊鏈(8255) C#(7972) AI(7469) 爪哇(7425) MySQL(7132) html(6777) 基礎類(6313) sql(6102) 熊猫(6058) PHP(5869) 数组(5741) R(5409) Linux(5327) 反应(5209) 腳本語言(PerlPython)(5129) 非技術區(4971) Android(4554) 数据框(4311) css(4259) 节点.js(4032) C語言(3288) json(3245) 列表(3129) 扑(3119) C++語言(3117) 安卓(2998) 打字稿(2995) VBA(2789) Java相關(2746) 疑難問題(2699) 细绳(2522) 單片機工控(2479) iOS(2429) ASP.NET(2402) MongoDB(2323) 麻木的(2285) 正则表达式(2254) 字典(2211) 循环(2198) 迅速(2185) 擅长(2169) 镖(2155) 功能(1967) .NET技术(1958) Web開發(1951) python-3.x(1918) HtmlCss(1915) 弹簧靴(1913) C++(1909) xml(1889) PostgreSQL(1872) .NETCore(1853) 谷歌表格(1846) Unity3D(1843) for循环(1842)

熱門瀏覽
  • 【C++】Microsoft C++、C 和匯編程式檔案

    ......

    uj5u.com 2020-09-10 00:57:23 more
  • 例外宣告

    相比于斷言適用于排除邏輯上不可能存在的狀態,例外通常是用于邏輯上可能發生的錯誤。 例外宣告 Item 1:當函式不可能拋出例外或不能接受拋出例外時,使用noexcept 理由 如果不打算拋出例外的話,程式就會認為無法處理這種錯誤,并且應當盡早終止,如此可以有效地阻止例外的傳播與擴散。 示例 //不可 ......

    uj5u.com 2020-09-10 00:57:27 more
  • Codeforces 1400E Clear the Multiset(貪心 + 分治)

    鏈接:https://codeforces.com/problemset/problem/1400/E 來源:Codeforces 思路:給你一個陣列,現在你可以進行兩種操作,操作1:將一段沒有 0 的區間進行減一的操作,操作2:將 i 位置上的元素歸零。最終問:將這個陣列的全部元素歸零后操作的最少 ......

    uj5u.com 2020-09-10 00:57:30 more
  • UVA11610 【Reverse Prime】

    本人看到此題沒有翻譯,就附帶了一個自己的翻譯版本 思考 這一題,它的第一個要求是找出所有 $7$ 位反向質數及其質因數的個數。 我們應該需要質數篩篩選1~$10^{7}$的所有數,這里就不慢慢介紹了。但是,重讀題,我們突然發現反向質數都是 $7$ 位,而將它反過來后的數字卻是 $6$ 位數,這就說明 ......

    uj5u.com 2020-09-10 00:57:36 more
  • 統計區間素數數量

    1 #pragma GCC optimize(2) 2 #include <bits/stdc++.h> 3 using namespace std; 4 bool isprime[1000000010]; 5 vector<int> prime; 6 inline int getlist(int ......

    uj5u.com 2020-09-10 00:57:47 more
  • C/C++編程筆記:C++中的 const 變數詳解,教你正確認識const用法

    1、C中的const 1、區域const變數存放在堆疊區中,會分配記憶體(也就是說可以通過地址間接修改變數的值)。測驗代碼如下: 運行結果: 2、全域const變數存放在只讀資料段(不能通過地址修改,會發生寫入錯誤), 默認為外部聯編,可以給其他源檔案使用(需要用extern關鍵字修飾) 運行結果: ......

    uj5u.com 2020-09-10 00:58:04 more
  • 【C++犯錯記錄】VS2019 MFC添加資源不懂如何修改資源宏ID

    1. 首先在資源視圖中,添加資源 2. 點擊新添加的資源,復制自動生成的ID 3. 在解決方案資源管理器中找到Resource.h檔案,編輯,使用整個專案搜索和替換的方式快速替換 宏宣告 4. Ctrl+Shift+F 全域搜索,點擊查找全部,然后逐個替換 5. 為什么使用搜索替換而不使用屬性視窗直 ......

    uj5u.com 2020-09-10 00:59:11 more
  • 【C++犯錯記錄】VS2019 MFC不懂的批量添加資源

    1. 打開資源頭檔案Resource.h,在其中預先定義好宏 ID(不清楚其實ID值應該設定多少,可以先新建一個相同的資源項,再在這個資源的ID值的基礎上遞增即可) 2. 在資源視圖中選中專案資源,按F7編輯資源檔案,按 ID 型別 相對路徑的形式添加 資源。(別忘了先把檔案拷貝到專案中的res檔案 ......

    uj5u.com 2020-09-10 01:00:19 more
  • C/C++編程筆記:關于C++的參考型別,專供新手入門使用

    今天要講的是C++中我最喜歡的一個用法——參考,也叫別名。 參考就是給一個變數名取一個變數名,方便我們間接地使用這個變數。我們可以給一個變數創建N個參考,這N + 1個變數共享了同一塊記憶體區域。(參考型別的變數會占用記憶體空間,占用的記憶體空間的大小和指標型別的大小是相同的。雖然參考是一個物件的別名,但 ......

    uj5u.com 2020-09-10 01:00:22 more
  • 【C/C++編程筆記】從頭開始學習C ++:初學者完整指南

    眾所周知,C ++的學習曲線陡峭,但是花時間學習這種語言將為您的職業帶來奇跡,并使您與其他開發人員區分開。您會更輕松地學習新語言,形成真正的解決問題的技能,并在編程的基礎上打下堅實的基礎。 C ++將幫助您養成良好的編程習慣(即清晰一致的編碼風格,在撰寫代碼時注釋代碼,并限制類內部的可見性),并且由 ......

    uj5u.com 2020-09-10 01:00:41 more
最新发布
  • Rust中的智能指標:Box<T> Rc<T> Arc<T> Cell<T> RefCell<T> Weak

    Rust中的智能指標是什么 智能指標(smart pointers)是一類資料結構,是擁有資料所有權和額外功能的指標。是指標的進一步發展 指標(pointer)是一個包含記憶體地址的變數的通用概念。這個地址參考,或 ” 指向”(points at)一些其 他資料 。參考以 & 符號為標志并借用了他們所 ......

    uj5u.com 2023-04-20 07:24:10 more
  • Java的值傳遞和參考傳遞

    值傳遞不會改變本身,參考傳遞(如果傳遞的值需要實體化到堆里)如果發生修改了會改變本身。 1.基本資料型別都是值傳遞 package com.example.basic; public class Test { public static void main(String[] args) { int ......

    uj5u.com 2023-04-20 07:24:04 more
  • [2]SpinalHDL教程——Scala簡單入門

    第一個 Scala 程式 shell里面輸入 $ scala scala> 1 + 1 res0: Int = 2 scala> println("Hello World!") Hello World! 檔案形式 object HelloWorld { /* 這是我的第一個 Scala 程式 * 以 ......

    uj5u.com 2023-04-20 07:23:58 more
  • 理解函式指標和回呼函式

    理解 函式指標 指向函式的指標。比如: 理解函式指標的偽代碼 void (*p)(int type, char *data); // 定義一個函式指標p void func(int type, char *data); // 宣告一個函式func p = func; // 將指標p指向函式func ......

    uj5u.com 2023-04-20 07:23:52 more
  • Django筆記二十五之資料庫函式之日期函式

    本文首發于公眾號:Hunter后端 原文鏈接:Django筆記二十五之資料庫函式之日期函式 日期函式主要介紹兩個大類,Extract() 和 Trunc() Extract() 函式作用是提取日期,比如我們可以提取一個日期欄位的年份,月份,日等資料 Trunc() 的作用則是截取,比如 2022-0 ......

    uj5u.com 2023-04-20 07:23:45 more
  • 一天吃透JVM面試八股文

    什么是JVM? JVM,全稱Java Virtual Machine(Java虛擬機),是通過在實際的計算機上仿真模擬各種計算機功能來實作的。由一套位元組碼指令集、一組暫存器、一個堆疊、一個垃圾回收堆和一個存盤方法域等組成。JVM屏蔽了與作業系統平臺相關的資訊,使得Java程式只需要生成在Java虛擬機 ......

    uj5u.com 2023-04-20 07:23:31 more
  • 使用Java接入小程式訂閱訊息!

    更新完微信服務號的模板訊息之后,我又趕緊把微信小程式的訂閱訊息給實作了!之前我一直以為微信小程式也是要企業才能申請,沒想到小程式個人就能申請。 訊息推送平臺🔥推送下發【郵件】【短信】【微信服務號】【微信小程式】【企業微信】【釘釘】等訊息型別。 https://gitee.com/zhongfuch ......

    uj5u.com 2023-04-20 07:22:59 more
  • java -- 緩沖流、轉換流、序列化流

    緩沖流 緩沖流, 也叫高效流, 按照資料型別分類: 位元組緩沖流:BufferedInputStream,BufferedOutputStream 字符緩沖流:BufferedReader,BufferedWriter 緩沖流的基本原理,是在創建流物件時,會創建一個內置的默認大小的緩沖區陣列,通過緩沖 ......

    uj5u.com 2023-04-20 07:22:49 more
  • Java-SpringBoot-Range請求頭設定實作視頻分段傳輸

    老實說,人太懶了,現在基本都不喜歡寫筆記了,但是網上有關Range請求頭的文章都太水了 下面是抄的一段StackOverflow的代碼...自己大修改過的,寫的注釋挺全的,應該直接看得懂,就不解釋了 寫的不好...只是希望能給視頻網站開發的新手一點點幫助吧. 業務場景:視頻分段傳輸、視頻多段傳輸(理 ......

    uj5u.com 2023-04-20 07:22:42 more
  • Windows 10開發教程_編程入門自學教程_菜鳥教程-免費教程分享

    教程簡介 Windows 10開發入門教程 - 從簡單的步驟了解Windows 10開發,從基本到高級概念,包括簡介,UWP,第一個應用程式,商店,XAML控制元件,資料系結,XAML性能,自適應設計,自適應UI,自適應代碼,檔案管理,SQLite資料庫,應用程式到應用程式通信,應用程式本地化,應用程式 ......

    uj5u.com 2023-04-20 07:22:35 more