【学习笔记】杜教筛

1.莫比乌斯求和

$$\begin{array}{rl}ϵ(n)& =\{\begin{array}{l}1,n=1\\ 0,n\ne 1\end{array}\\ I(n)& =1\\ ϵ(n)& =\sum _{i\mid n}\mu (i)\times I(\frac{n}{i})\\ \sum _{i=1}^nϵ(i)& =\sum _{i=1}^n\sum _{d\mid i}\mu (d)\times I(\frac{i}{d})\\ & =\sum _{i=1}^{n}I(i)\times Q(\lfloor \frac{n}{i}\rfloor )\\ I(1)\times Q(n)& =\sum _{i=1}^nϵ(i)-\sum _{i=2}^{n}I(i)\times Q(\lfloor \frac{n}{i}\rfloor )\end{array}$$

#include <map>
#include <cmath>
#include <queue>
#include <vector>
#include <cstdio>
#include <cstdlib> 
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define fi first
#define se second
#define db double
#define LL long long
#define ULL unsigned long long
#define PII pair <int, int>
#define MP(x,y) make_pair (x, y)
#define rep(i,j,k) for (int i = (j); i <= (k); ++i)
#define per(i,j,k) for (int i = (j); i >= (k); --i)

template <typename T> T Max (T x, T y) { return x > y ? x : y; }
template <typename T> T Min (T x, T y) { return x < y ? x : y; }
template <typename T> T Abs (T x) { return x > 0 ? x : -x; }
template <typename T>
void read (T &x) {
    x = 0; T f = 1;
    char ch = getchar ();
    while (ch < '0' || ch > '9') {
        if (ch == '-') f = -1;
        ch = getchar ();
    }
    while (ch >= '0' && ch <= '9') {
        x = (x << 3) + (x << 1) + ch - '0';
        ch = getchar ();
    }
    x *= f;
}
template <typename T, typename... Args>
void read (T &x, Args&... args) {
    read (x); read (args...);
}
char For_Print[25];
template <typename T>
void write (T x) {
    if (x == 0) { putchar ('0'); return; }
    if (x < 0) { putchar ('-'); x = -x; }
    int poi = 0;
    while (x) {
        For_Print[++poi] = x % 10 + '0';
        x /= 10;
    }
    while (poi) putchar (For_Print[poi--]);
}
template <typename T>
void print (T x, char ch) {
    write (x); putchar (ch);
}

const LL Mod = 1e9 + 7;

LL square (LL x) { return (x * x) % Mod; }
void DEL (LL &x, LL y) { ((x -= y) < 0) && (x += Mod); }
void ADD (LL &x, LL y) { ((x += y) >= Mod) && (x -= Mod); }

const int Maxn = 1e6;

LL t, l, r;

bool vis[Maxn + 5];
int cnt, primes[Maxn + 5];
int mu[Maxn + 5], pre[Maxn + 5];
void Euler () {
	mu[1] = 1;
	for (int i = 2; i <= Maxn; i++) {
		if (vis[i] == 0) {
			primes[++cnt] = i;
			mu[i] = -1;
		}
		rep (j, 1, cnt) {
			if (primes[j] > Maxn / i) break;
			vis[i * primes[j]] = 1;
			if (i % primes[j] == 0) {
				mu[i * primes[j]] = 0;
				break;
			}
			mu[i * primes[j]] = -mu[i];
		}
	}
	rep (i, 1, Maxn) {
		pre[i] = pre[i - 1] + mu[i];
	}
}
map <LL, int> dp;
LL Calc (LL n) {
	if (n <= Maxn) return pre[n];
	if (dp.find (n) != dp.end ()) return dp[n];
	LL l = 2, r, res = 1;
	while (l <= n) {
		r = (n / (n / l));
		res -= (r - l + 1) * Calc (n / l);
		l = r + 1;
	}
	return dp[n] = res;
}

int main () {
//	freopen ("D:\\lihan\\1.in", "r", stdin);
//	freopen ("D:\\lihan\\1.out", "w", stdout);

	Euler ();

	read (l, r);
	cout << Calc (r) - Calc (l - 1) << endl;
    return 0;
}

---

2.欧拉函数求和

$$\begin{array}{rl}id(n)& =\sum _{d\mid n}\phi (d)\\ \sum _{i=1}^{n}id(i)& =\sum _{i=1}^n\sum _{d\mid i}\phi (d)\\ & =\sum _{i=1}^n\sum _{d\mid i}\phi (d)I(\frac{n}{d})\\ & =\sum _{i=1}^{n}I(i)Q(\lfloor \frac{n}{i}\rfloor )\\ I(1)Q(n)& =\sum _{i=1}^{n}id(i)-\sum _{i=2}^{n}I(i)Q(\lfloor \frac{n}{i}\rfloor )\end{array}$$

#include <map>
#include <cmath>
#include <queue>
#include <vector>
#include <cstdio>
#include <cstdlib> 
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define fi first
#define se second
#define db double
#define LL long long
#define ULL unsigned long long
#define PII pair <int, int>
#define MP(x,y) make_pair (x, y)
#define rep(i,j,k) for (int i = (j); i <= (k); ++i)
#define per(i,j,k) for (int i = (j); i >= (k); --i)

template <typename T> T Max (T x, T y) { return x > y ? x : y; }
template <typename T> T Min (T x, T y) { return x < y ? x : y; }
template <typename T> T Abs (T x) { return x > 0 ? x : -x; }
template <typename T>
void read (T &x) {
    x = 0; T f = 1;
    char ch = getchar ();
    while (ch < '0' || ch > '9') {
        if (ch == '-') f = -1;
        ch = getchar ();
    }
    while (ch >= '0' && ch <= '9') {
        x = (x << 3) + (x << 1) + ch - '0';
        ch = getchar ();
    }
    x *= f;
}
template <typename T, typename... Args>
void read (T &x, Args&... args) {
    read (x); read (args...);
}
char For_Print[25];
template <typename T>
void write (T x) {
    if (x == 0) { putchar ('0'); return; }
    if (x < 0) { putchar ('-'); x = -x; }
    int poi = 0;
    while (x) {
        For_Print[++poi] = x % 10 + '0';
        x /= 10;
    }
    while (poi) putchar (For_Print[poi--]);
}
template <typename T>
void print (T x, char ch) {
    write (x); putchar (ch);
}

const int Maxn = 1e6;
const LL Mod = 1000000007;

LL m;

bool vis[Maxn + 5];
int cnt, primes[Maxn + 5];
LL varphi[Maxn + 5], pre[Maxn + 5];
void Euler () {
	varphi[1] = 1;
	for (int i = 2; i <= Maxn; i++) {
		if (vis[i] == 0) {
			primes[++cnt] = i;
			varphi[i] = i - 1;
		}
		rep (j, 1, cnt) {
			if (primes[j] > Maxn / i) break;
			vis[i * primes[j]] = 1;
			if (i % primes[j] == 0) {
				varphi[i * primes[j]] = varphi[i] * primes[j];
				break;
			}
			varphi[i * primes[j]] = varphi[i] * (primes[j] - 1);
		}
	}
	rep (i, 1, Maxn) {
		pre[i] = pre[i - 1] + varphi[i];
		pre[i] %= Mod;
	}
}
map <LL, LL> dp;
LL Calc (LL n) {
	if (n <= Maxn) return pre[n];
	if (dp.find (n) != dp.end ()) return dp[n];
	LL l = 2, r, res = ((n % Mod) * ((n + 1) % Mod) / 2) % Mod;
	while (l <= n) {
		r = (n / (n / l));
		res = (res - Calc (n / l) * ((r - l + 1) % Mod) % Mod + Mod) % Mod;
		l = r + 1;
	}
	return dp[n] = res;
}

int main () {
//	freopen ("D:\\lihan\\1.in", "r", stdin);
//	freopen ("D:\\lihan\\1.out", "w", stdout);

	Euler ();

	read (m);
	write (Calc (m));
    return 0;
}

---

3.欧拉函数乘上 $id$

$$\begin{array}{rl}Q(i)& =\sum _{i=1}^{n}i\times \phi (i)\\ (Q\ast id)(n)& =\sum _{d\mid n}i\times \phi (i)\times id(\frac{n}{i})\\ & =\sum _{d\mid n}\phi (i)\times n\\ & =n\times n\\ \sum _{i=1}^n{i}^2& =\sum _{i=1}^n\sum _{d\mid i}d\times \phi (d)\times id(\frac{i}{d})\\ & =\sum _{i=1}^{n}id(i)\times Q(\lfloor \frac{n}{i}\rfloor )\\ Q(n)& =\sum _{i=1}^n{i}^2-\sum _{i=2}^{n}id(i)\times Q(\lfloor \frac{n}{i}\rfloor )\end{array}$$

---

4.数论函数求和

$$\begin{array}{rl}\sum _{d\mid n}f(d)& =n^2-3n+2\\ Q(n)& =\sum _{i=1}^{n}f(i)\\ \sum _{i=1}^n(i^2-3i+2)& =\sum _{i=1}^n\sum _{d\mid i}f(d)I(\frac{n}{d})\\ & =\sum _{i=1}^{n}I(i)\times Q(\lfloor \frac{n}{i}\rfloor )\\ Q(n)& =\sum _{i=1}(i^2-3i+2)-\sum _{i=2}^{n}I(i)\times Q(\lfloor \frac{n}{i}\rfloor )\end{array}$$

---

5.平均最小公倍数

引理 $1$:

${\sum }_{i=1}^n\frac{lcm(i,n)}{n}=\frac{1+{\sum }_{d\mid n}\phi (d)d}{2}$

证明：

$$\begin{array}{rl}\sum _{i=1}^n\frac{lcm(i,n)}{n}& =\sum _{g=1}^n\sum _{i=1}^n\frac{i\times n}{g\times n}[g=gcd(i,n)]\\ & =\sum _{g=1}^n\sum _{i=1}^n\frac{i}{g}\cdot [g=gcd(i,n)]\\ & =\sum _{g\mid n}\sum _{i=1}^n\frac{i}{g}\cdot [g=gcd(i,n)]\\ & =\sum _{g\mid n}\sum _{i=1}^n\frac{i}{g}\cdot [g=gcd(i,n)]\end{array}$$

---

$$\begin{array}{rl}Q(n)& =\sum _{i=1}^n\sum _{j=1}^i\frac{lcm(i,j)}{i}\\ Q(n)& =\sum _{i=1}^n\frac{1+{\sum }_{d\mid i}\phi (d)d}{2}(引理)\\ Q(n)& =\frac{1}{2}(\sum _{i=1}^n\sum _{d\mid i}\phi (d)d)+\frac{n}{2}\\ 即求& \sum _{i=1}^n\sum _{d\mid i}\phi (d)d\\ \sum _{i=1}^n\sum _{d\mid i}\phi (d)d& =\sum _{i=1}^n\phi (i)\times i\times \lfloor \frac{n}{i}\rfloor \\ 数论& 分块，即求(\sum _{i=l}^r\phi (i)\times i)\times \lfloor \frac{n}{l}\rfloor （用3.）\end{array}$$

外面套了一个数论分块，似乎时间复杂度是 $O(n^{\frac{1}{2}}\cdot n^{\frac{1}{2}})=O(n)$

但是询问的 $m$ 都是满足 $\exists i,\lfloor \frac{n}{i}\rfloor$，所以询问不会超过 $\sqrt{(}n)$ 个。

所以时间复杂度还是没变。

#include <map>
#include <cmath>
#include <queue>
#include <vector>
#include <cstdio>
#include <cstdlib> 
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define fi first
#define se second
#define db double
#define LL long long
#define ULL unsigned long long
#define PII pair <int, int>
#define MP(x,y) make_pair (x, y)
#define rep(i,j,k) for (int i = (j); i <= (k); ++i)
#define per(i,j,k) for (int i = (j); i >= (k); --i)

template <typename T> T Max (T x, T y) { return x > y ? x : y; }
template <typename T> T Min (T x, T y) { return x < y ? x : y; }
template <typename T> T Abs (T x) { return x > 0 ? x : -x; }
template <typename T>
void read (T &x) {
    x = 0; T f = 1;
    char ch = getchar ();
    while (ch < '0' || ch > '9') {
        if (ch == '-') f = -1;
        ch = getchar ();
    }
    while (ch >= '0' && ch <= '9') {
        x = (x << 3) + (x << 1) + ch - '0';
        ch = getchar ();
    }
    x *= f;
}
template <typename T, typename... Args>
void read (T &x, Args&... args) {
    read (x); read (args...);
}
char For_Print[25];
template <typename T>
void write (T x) {
    if (x == 0) { putchar ('0'); return; }
    if (x < 0) { putchar ('-'); x = -x; }
    int poi = 0;
    while (x) {
        For_Print[++poi] = x % 10 + '0';
        x /= 10;
    }
    while (poi) putchar (For_Print[poi--]);
}
template <typename T>
void print (T x, char ch) {
    write (x); putchar (ch);
}

const int Maxn = 1e6;
const LL Mod = 1000000007;

LL md (LL x, LL M) {
	return (x % M + M) % M;
}

LL inv_2, inv_6;

LL quick_pow (LL x, LL y) {
	LL res = 1;
	while (y) {
		if (y & 1) res = md (res * x, Mod);
		x = md (x * x, Mod); y >>= 1;
	}
	return res;
}
bool vis[Maxn + 5];
int cnt, primes[Maxn + 5];
LL phi[Maxn + 5], pre[Maxn + 5];
void Euler () {
	phi[1] = 1;
	for (int i = 2; i <= Maxn; i++) {
		if (vis[i] == 0) {
			primes[++cnt] = i;
			phi[i] = i - 1;
		}
		rep (j, 1, cnt) {
			if (primes[j] > Maxn / i) break;
			vis[i * primes[j]] = 1;
			if (i % primes[j] == 0) {
				phi[i * primes[j]] = phi[i] * primes[j];
				break;
			}
			phi[i * primes[j]] = phi[i] * (primes[j] - 1);
		}
	}
	rep (i, 1, Maxn) {
		pre[i] = pre[i - 1] + md (phi[i] * i, Mod);
		pre[i] %= Mod;
	}
}
map <LL, LL> dp;
LL Calc (LL n) {
	if (n <= Maxn) return pre[n];
	if (dp.find (n) != dp.end ()) return dp[n];
	LL l = 2, r, res = md (md (n, Mod) * md (md (n + 1, Mod) * md (md (2 * n + 1, Mod) * inv_6, Mod), Mod), Mod);
	while (l <= n) {
		r = (n / (n / l));
		res = md (res - md (Calc (n / l) * md (md (l + r, Mod) * md (md (r - l + 1, Mod) * inv_2, Mod), Mod), Mod), Mod);
		l = r + 1;
	}
	return dp[n] = res;
}
LL Query (LL n) {
	LL l = 1, r, res = 0;
	while (l <= n) {
		r = n / (n / l);
		res = md (res + md (Calc (r) - Calc (l - 1), Mod) * (n / l), Mod);
		l = r + 1;
	}
	res = md (res * inv_2, Mod);
	res = md (res + md (n * inv_2, Mod), Mod);
	return res;
}

int main () {
//	freopen ("D:\\lihan\\1.in", "r", stdin);
//	freopen ("D:\\lihan\\1.out", "w", stdout);

	Euler ();
	
	inv_2 = quick_pow (2, Mod - 2);
	inv_6 = quick_pow (6, Mod - 2);
	LL l, r; cin >> l >> r;
	cout << md (Query (r) - Query (l - 1), Mod);
    return 0;
}

---

6.Lucas的数论

引理：

$$\sum _{d\mid ij}1(d)=\sum _{x\mid i}\sum _{y\mid j}[gcd(x,y)=1]$$

算二次证明。

相当于我们枚举的因数为： $\frac{i}{x}\cdot y$。

不漏：

每个 $ij$ 的因数都应该是枚举得到的，不用哆嗦。

不重：

相当于我们的 $x$ 会从 $i$ 处刮掉一些因子，而 $y$ 会在 $j$ 中取一些因子，$gcd(x,y)=1$ 就防止了从 $i$ 处刮掉的因子和从 $j$ 处取得的因子出现重复，就不会算重了。

---

$$\begin{array}{rl}\sum _{i=1}^n\sum _{j=1}^{n}f(ij)& =\sum _{i=1}^n\sum _{j=1}^n\sum _{x\mid i}\sum _{y\mid j}[gcd(x,y)=1]\\ & =\sum _{i=1}^n\sum _{j=1}^n\sum _{x\mid i}\sum _{y\mid j}\sum _{t\mid gcd(x,y)}\mu (t)\\ & =\sum _{x=1}^n\sum _{y=1}^n\sum _{x\mid i}^n\sum _{y\mid j}^n\sum _{t\mid gcd(x,y)}\mu (t)\\ & =\sum _{x=1}^n\sum _{y=1}^n\lfloor \frac{n}{x}\rfloor \lfloor \frac{n}{y}\rfloor \sum _{t\mid gcd(x,y)}\mu (t)\\ & =\sum _{t=1}^n\mu (t)\sum _{t\mid x}\lfloor \frac{n}{x}\rfloor \sum _{t\mid y}\lfloor \frac{n}{y}\rfloor \\ & =\sum _{t=1}^n\mu (t)\sum _{x=1}^{\lfloor \frac{n}{t}\rfloor }\lfloor \frac{\lfloor \frac{n}{t}\rfloor }{x}\rfloor \sum _{y=1}^{\lfloor \frac{n}{t}\rfloor }\lfloor \frac{\lfloor \frac{n}{t}\rfloor }{y}\rfloor \\ 令f(n)& =\sum _{i=1}^n\lfloor \frac{n}{i}\rfloor \\ 原式& =\sum _{t=1}^n\mu (t)f(\lfloor \frac{n}{t}\rfloor )f(\lfloor \frac{n}{t}\rfloor )\end{array}$$

在外面一个数论分块，$f(\lfloor \frac{n}{t}\rfloor )$ 可以直接暴力求（根据杜教筛的时间复杂度理论），然后 $\mu (t)$ 求和直接用杜教筛就可以啦。

#include <map>
#include <cmath>
#include <queue>
#include <vector>
#include <cstdio>
#include <cstdlib> 
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define fi first
#define se second
#define db double
#define LL long long
#define ULL unsigned long long
#define PII pair <int, int>
#define MP(x,y) make_pair (x, y)
#define rep(i,j,k) for (int i = (j); i <= (k); ++i)
#define per(i,j,k) for (int i = (j); i >= (k); --i)

template <typename T> T Max (T x, T y) { return x > y ? x : y; }
template <typename T> T Min (T x, T y) { return x < y ? x : y; }
template <typename T> T Abs (T x) { return x > 0 ? x : -x; }
template <typename T>
void read (T &x) {
    x = 0; T f = 1;
    char ch = getchar ();
    while (ch < '0' || ch > '9') {
        if (ch == '-') f = -1;
        ch = getchar ();
    }
    while (ch >= '0' && ch <= '9') {
        x = (x << 3) + (x << 1) + ch - '0';
        ch = getchar ();
    }
    x *= f;
}
template <typename T, typename... Args>
void read (T &x, Args&... args) {
    read (x); read (args...);
}
char For_Print[25];
template <typename T>
void write (T x) {
    if (x == 0) { putchar ('0'); return; }
    if (x < 0) { putchar ('-'); x = -x; }
    int poi = 0;
    while (x) {
        For_Print[++poi] = x % 10 + '0';
        x /= 10;
    }
    while (poi) putchar (For_Print[poi--]);
}
template <typename T>
void print (T x, char ch) {
    write (x); putchar (ch);
}

const int Maxn = 1e6;
const LL Mod = 1000000007;

LL md ( LL x ) {
    return ( x % Mod + Mod ) % Mod;
}

LL t, l, r;

bool vis[Maxn + 5];
int cnt, primes[Maxn + 5];
LL mu[Maxn + 5], pre[Maxn + 5];
void Euler () {
	mu[1] = 1;
	for (int i = 2; i <= Maxn; i++) {
		if (vis[i] == 0) {
			primes[++cnt] = i;
			mu[i] = -1;
		}
		rep (j, 1, cnt) {
			if (primes[j] > Maxn / i) break;
			vis[i * primes[j]] = 1;
			if (i % primes[j] == 0) {
				mu[i * primes[j]] = 0;
				break;
			}
			mu[i * primes[j]] = -mu[i];
		}
	}
	rep (i, 1, Maxn) {
		pre[i] = pre[i - 1] + mu[i];
	}
}
map <LL, LL> dp;
LL Calc (LL n) {
	if (n <= Maxn) return pre[n];
	if (dp.find (n) != dp.end ()) return dp[n];
	LL l = 2, r, res = 1;
	while (l <= n) {
		r = (n / (n / l));
		res = md (res - md ( md (r - l + 1) * md ( Calc ( n / l ) ) ) );
		l = r + 1;
	}
	return dp[n] = res;
}

LL F ( LL x ) {
    LL l = 1, r, res = 0;
    while ( l <= x ) {
        r = x / ( x / l );
        res = md ( res + md ( r - l + 1 ) * md ( x / l ) );
        l = r + 1;
    }
    return res;
}

LL n;

int main () {
	// freopen ( "C:\\Users\\cqbz\\Desktop\\lihan\\1.in", "r", stdin );
	// freopen ( "C:\\Users\\cqbz\\Desktop\\lihan\\1.out", "w", stdout );

	Euler ();

    read ( n );
    Calc ( n );

    LL l = 1, r, res = 0;
    while ( l <= n ) {
        r = n / ( n / l );
        LL tmp = F ( n / l );
        res = md ( res + md ( md ( md ( Calc ( r ) - Calc ( l - 1 ) ) * tmp ) * tmp ) );
        l = r + 1;
    }
    write ( res );
    return 0;
}

7. 约数之和

$$\begin{array}{rl}\sum _{i=1}^n\sum _{j=1}^{n}f(ij)& =\sum _{i=1}^n\sum _{j=1}^n\sum _{x\mid i}\sum _{y\mid j}[gcd(x,y)=1]\cdot \frac{i}{x}\cdot y\\ & =\sum _{i=1}^n\sum _{j=1}^n\sum _{x\mid i}\sum _{y\mid j}\sum _{t\mid gcd(x,y)}\mu (t)\cdot \frac{i}{x}\cdot y\\ & =\sum _{x=1}^n\sum _{y=1}^n\sum _{x\mid i}\sum _{y\mid j}\sum _{t\mid x,t\mid y}\mu (t)\cdot (i/x)\cdot y\\ & =\sum _{x=1}^n\sum _{y=1}^n\sum _{x\mid i}\sum _{y\mid j}\sum _{t\mid x,t\mid y}\mu (t)\cdot (i/x)\cdot y\\ & =\sum _{x=1}^n\sum _{y=1}^{n}y\sum _{x\mid i}(i/x)\sum _{y\mid j}\sum _{t\mid x,t\mid y}\mu (t)\\ & =\sum _t\mu (t)\sum _{t\mid x}^n\sum _{t\mid y}^{n}y\sum _{x\mid i}^n(i/x)\sum _{y\mid j}^{n}1\\ & =\sum _t\mu (t)\sum _x^{\lfloor \frac{n}{t}\rfloor }\sum _y^{\lfloor \frac{n}{t}\rfloor }yt\sum _{xt\mid i}(i/xt)\sum _{yt\mid j}1\\ & =\sum _t\mu (t)\sum _x^{\lfloor \frac{n}{t}\rfloor }\sum _{xt\mid i}(i/xt)\sum _y^{\lfloor \frac{n}{t}\rfloor }yt\sum _{yt\mid j}1\\ & =\sum _t\mu (t)\sum _x^{\lfloor \frac{n}{t}\rfloor }(1+\lfloor \frac{\lfloor \frac{n}{t}\rfloor }{x}\rfloor )\ast \lfloor \frac{\lfloor \frac{n}{t}\rfloor }{x}\rfloor /2\cdot t\sum _y^{\lfloor \frac{n}{t}\rfloor }y\lfloor \frac{\lfloor \frac{n}{t}\rfloor }{y}\rfloor \\ g(n)& =\sum _{x=1}^n(1+\lfloor \frac{n}{x}\rfloor )\lfloor \frac{n}{x}\rfloor /2\\ f(n)& =\sum _{x=1}^n\lfloor \frac{n}{x}\rfloor x\\ & =\sum _t\mu (t)\cdot t\cdot g(\lfloor \frac{n}{t}\rfloor )\cdot f(\lfloor \frac{n}{t}\rfloor )\end{array}$$

这个式子不太好做 😓。

于是大部分做法都是把 $i,j$ 用 $x,y$ 代掉。

$$\begin{array}{rl}\sum _{i=1}^n\sum _{j=1}^{n}f(ij)& =\sum _{i=1}^n\sum _{j=1}^n\sum _{x\mid i}\sum _{y\mid j}[gcd(x,y)=1]\cdot \frac{i}{x}\cdot y\\ & =\sum _{i=1}^n\sum _{j=1}^n\sum _{x\mid i}\sum _{y\mid j}\sum _{t\mid gcd(x,y)}\mu (t)\cdot \frac{i}{x}\cdot y\\ & =\sum _{x=1}^n\sum _{y=1}^n\sum _{x\mid i}\sum _{y\mid j}\sum _{t\mid x,t\mid y}\mu (t)\cdot (i/x)\cdot y\\ & =\sum _{x=1}^n\sum _{y=1}^n\sum _{x\mid i}\sum _{y\mid j}\sum _{t\mid x,t\mid y}\mu (t)\cdot (i/x)\cdot y\\ & =\sum _{x=1}^n\sum _{y=1}^{n}y\sum _{x\mid i}(i/x)\sum _{y\mid j}\sum _{t\mid x,t\mid y}\mu (t)\\ & =\sum _t\mu (t)\sum _{t\mid x}^n\sum _{t\mid y}^{n}y\sum _{x\mid i}^n(i/x)\sum _{y\mid j}^{n}1\\ & =\sum _t\mu (t)\sum _{t\mid x}^n\sum _{t\mid y}^{n}y\sum _i^{\lfloor \frac{n}{i}\rfloor }i\lfloor \frac{n}{y}\rfloor \end{array}$$