DZY Loves Math 系列详细题解

BZOJ 3309: DZY Loves Math I

题意

(f(n)) 为 (n) 幂指数的最大值。

[sum_{i = 1}^{a} sum_{j = 1}^{b} f(gcd(i, j)) ]

(Tle 10000, 1 le a,b le 10^7)

题解

[egin{aligned} ans &= sum_{i = 1}^{a} sum_{j = 1}^{b} f(gcd(i, j)) \ &= sum_{d = 1}^{min(a, b)} f(d) sum_{x = 1}^{min(lfloor frac a d floor, lfloor frac b d floor)} mu(x) lfloor frac a {dx} floor lfloor frac b {dx} floor end{aligned} ]

令 (dx = T) 那么有

[egin{aligned} ans &= sum_{T = 1}^{min(a, b)} lfloor frac a {T} floor lfloor frac b {T} floor (sum_{d | T} f(d) mu(frac T d)) end{aligned} ]

我们线性筛出 (f * mu) 就行了，每次整除分块回答，复杂度是 (mathcal O(n + T sqrt n)) 的。

什么不会线性筛？那么埃氏筛卡常吧。

BZOJ 3462: DZY Loves Math II

题意

(2le S le 2*10^6,1 le n le 10^{18},1 le q le 10^5)

题解

好神啊。。

我们考虑把相同的 (p_i) 合并，那么就变成

[egin{aligned} sum_{i le k} p_i c_i &= n &(c_i ge 1) \ prod_{i le k} p_i &= S end{aligned} ]

求第一个方程 (c_i) 解的方案数。那么对于 (S) 的唯一分解来说，不存在一个质因子的次数 (>1) 。

对于 (c_i ge 1) 的限制我们可以把 (n - sum_i p_i) 来去掉，那么现在我们只需要求一个完全背包的方案数。

对于一个因子 (p_i) 选 (c_i) 个的体积是 (p_i c_i) ，其实就等价于 (xS + yp_i) 这个方案。

什么意思呢？我们把 (p_ic_i) 对于 (S) 取模后就得到了 (yp_i) ，也就意味着 (yp_i < S) 。

那么我们要求的其实就是

[egin{aligned} sum_{i le k} (xS + yp_i) &= n &(c_i ge 1, yp_i < S) \ end{aligned} ]

的方案数。

考虑移项，那么就变成

[egin{aligned} sum_{i le k} yp_i &= n - sum_{i le k} xS &(c_i ge 1, y < frac S {p_i}) \ end{aligned} ]

不难发现对于前者的大小不会超过 (kS) ，可以直接暴力做多重背包就行了。

那么对于后面的方案如何算呢？不难发现其实就是求 (sum_{i le k} x = displaystyle frac{n - sum_{i le k}y p_i}{S}) 的方案数，直接隔板一下就行了。

复杂度其实是 (mathcal O(k^2S + qk^2)) 的，跑的挺快。

总结

对于一类计数问题，可以考虑如何把原来模型等价替代，常常能优化复杂度，或减少难度。

代码

注意多重背包那里可以减掉不合法的优化一下复杂度。

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl

using namespace std;

typedef long long ll;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

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

void File() {
#ifdef zjp_shadow
	freopen ("3462.in", "r", stdin);
	freopen ("3462.out", "w", stdout);
#endif
}

const int N = 2e6 + 1e3, Mod = 1e9 + 7;

inline void add(int &a, int b) {
	if ((a += b) >= Mod) a -= Mod;
}

inline int fpm(int x, int power) {
	int res = 1;
	for (; power; power >>= 1, x = 1ll * x * x % Mod)
		if (power & 1) res = 1ll * res * x % Mod;
	return res;
}

int fac[N], ifac[N];
void Fac_Init(int maxn) {
	fac[0] = ifac[0] = 1;
	For (i, 1, maxn) fac[i] = 1ll * fac[i - 1] * i % Mod;
	ifac[maxn] = fpm(fac[maxn], Mod - 2);
	Fordown (i, maxn - 1, 1) ifac[i] = ifac[i + 1] * (i + 1ll) % Mod;
}

inline int comb(ll n, int m) {
	if (n < 0 || m < 0 || n < m) return 0;
	int res = 1;
	for (ll i = n; i >= n - m + 1; -- i)
		res = 1ll * (i % Mod) * res % Mod;
	return 1ll * res * ifac[m] % Mod;
}

int S, q, p[20], ptot, f[2][N * 7];

int main () {

	File();

	S = read(); int tmp = S; q = read();

	for (int i = 2; i * i <= tmp; ++ i)
		while (!(tmp % i)) tmp /= i, p[++ ptot] = i;
	if (tmp > 1) p[++ ptot] = tmp;
	sort(p + 1, p + ptot + 1);
	ptot = unique(p + 1, p + ptot + 1) - p - 1;
	if (p[ptot + 1]) {
		while (q --) puts("0"); return 0;
	}

	int cur = 0; f[0][0] = 1;
	For (i, 1, ptot) {
		cur ^= 1; Cpy(f[cur], f[cur ^ 1]);
		For (j, 1, i * S) {
			if (j >= p[i]) add(f[cur][j], f[cur][j - p[i]]);
			if (j >= S) add(f[cur][j], Mod - f[cur ^ 1][j - S]);
		}
	}

	int sum = accumulate(p + 1, p + ptot + 1, 0); Fac_Init(ptot);
	while (q --) {
		ll n = read() - sum, ans = 0;
		if (n < 0) {
			puts("0"); continue;
		}
		for (int i = n % S; i <= ptot * S; i += S)
			ans = (ans + 1ll * f[cur][i] * comb((n - i) / S + ptot - 1, ptot - 1)) % Mod;
		printf ("%lld
", ans);
	}

	return 0;

}

BZOJ 3481: DZY Loves Math III

题意

给定整数 (P, Q) 求满足方程 (xy equiv Q pmod P) 的整数解 ((x, y)) 的数量，满足 (0 le x, y < P) 。

其中 (P = prod_{i le N} P_i, Q = prod_{i le N} Q_i) 。

(1 le N le 10, 0 le Q_i le 10^{18}, 1 le P_i le 10^{18},P ge 2)

题解

考虑枚举一个 (x) ，那么我们就是求对于 (yx + kP = Q) 的整数解个数。

利用在扩欧学习的知识，只有在 (gcd(x, P) | Q) 的时候才会有解，且在这个模意义下有 (gcd(x, P)) 个解。

那么答案其实就是

[sum_{x < P} gcd(x, P) [gcd(x, P) | Q] ]

考虑枚举 (gcd(P, Q) = T) 那么答案其实就是

[f(T) = sum_{d | T} d cdot varphi(frac P d) ]

如果你足够厉害，就能发现 (f(T)) 其实是个积性函数。

其实是因为积性函数和完全积性函数的狄利克雷卷积是个积性函数。。。（把 (P) 那里提出一个因子变成 (T) 就行了）

那么我们就能对于每个质因子 (p) 单独考虑啦，设它在 (P) 中有 (q) 个，(T) 中有 (q') 个（显然有 (q' le q) ）。

然后贡献其实就是

[egin{aligned} ans &= sum_{i = 0}^{q'} p^i cdot varphi(p^{q - i})\ &= sum_{i = 0}^{q'} p^i cdot (p - 1) p^{q - i - 1}\ &= q' (p - 1) p^{q - 1} end{aligned} ]

但是这个会在 (q' = q) 的时候出现问题，因为 (varphi(1) = 1) 那么我们就会把个 (p^q) 算成 ((p - 1)p^{q - 1}) ，加回去即可。

总结

(yx + kP = Q) 有在 (gcd(x, P) | Q) 的时候才会有解，且在这个模 (P) 意义下有 (gcd(x, P)) 个解。

注意观察积性函数的性质，用质数去算答案大大降低复杂度。

代码

要个 ( ext{Pollard Rho}) 分解质因数qwq

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
#define pb push_back

using namespace std;

typedef long long ll;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

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

void File() {
#ifdef zjp_shadow
	freopen ("3481.in", "r", stdin);
	freopen ("3481.out", "w", stdout);
#endif
}

bool flag;

namespace Pollard_Rho {

	inline ll rnd() { return abs((rand() << 30ll) | rand()); }
	inline ll randint(ll l, ll r) { return rnd() % (r - l + 1) + l; }

	const double eps = 1e-3;
	inline ll mul(ll a, ll b, ll Mod) {
		ll tmp = (a * b - (ll) ((long double) a / Mod * b + eps) * Mod);
		return tmp < 0 ? tmp + Mod : tmp;
	}

	inline ll fpm(ll x, ll power, ll Mod) {
		ll res = 1;
		for(; power; power >>= 1, x = mul(x, x, Mod))
			if (power & 1) res = mul(res, x, Mod);
		return res;
	}

	const int times = 8;


	inline bool Miller_Rabin(ll p) {
		if (p <= 2) return p == 2;
		if (!(p & 1)) return false;
		ll u = p - 1; int power = 0;
		for (; !(u & 1); u >>= 1) ++ power;
		For (i, 1, times) {
			ll a = randint(2, p - 1), x = fpm(a, u, p), y;
			for (int i = 1; i <= power; ++ i, x = y) {
				if ((y = mul(x, x, p)) == 1 && x != 1 && x != p - 1) return false;
			}
			if (x != 1) return false;
		}
		return true;
	}

	ll c, Mod;

	ll f(ll x) { return (mul(x, x, Mod) + c) % Mod; }

	ll find(ll x) {
		if (!(x & 1)) return 2; Mod = x;
		ll a = randint(2, x - 1), b = a; c = randint(2, x - 1);
		do {
			a = f(a); b = f(f(b));
			ll p = __gcd(abs(a - b), x); 
			if (p > 1) return p;
		} while (b != a);
		return find(x);
	}

	void ReSolve(ll x, map<ll, int> &factor) {
		if (x <= 1) { flag |= !x; return; }
		if (Miller_Rabin(x)) { ++ factor[x]; return; }
		ll fac = find(x); ReSolve(fac, factor); ReSolve(x / fac, factor);
	}

}

map<ll, int> P, Q, T; typedef map<ll, int> :: iterator iter;

const int Mod = 1e9 + 7;

int main() {

	File();

	int n = read();

	srand(998244353);

	using Pollard_Rho :: ReSolve;
	using Pollard_Rho :: fpm;

	For (i, 1, n) ReSolve(read(), P);
	For (i, 1, n) ReSolve(read(), Q);
	if (flag) Q = P;

	for (iter it = Q.begin(); it != Q.end(); ++ it) {
		ll p = it -> first; T[p] = min(P[p], Q[p]);
	}

	int ans = 1;
	for (iter it = P.begin(); it != P.end(); ++ it) 
		if (it -> second) {
			ll p = it -> first, res;
			res = (p - 1) % Mod * (T[p] + 1) % Mod * fpm(p, P[p] - 1, Mod) % Mod;
			if (T[p] == P[p])
				res = (res - (p - 1) % Mod * fpm(p, P[p] - 1, Mod) + fpm(p, P[p], Mod)) % Mod;
			ans = 1ll * ans * (res + Mod) % Mod;
		}
	printf ("%d
", ans);

	return 0;

}

BZOJ 3512: DZY Loves Math IV

题意

给定 (n, m) 求

[sum_{i = 1}^{n} sum_{j = 1}^{m} varphi(ij) pmod {10^9+7} ]

(1 le n le 10^5, 1 le m le 10^9)

题解

神题啊！

为啥 (n, m) 不是同级的？因为是让我们暴力枚举 (n) 。

那么我们求得就是

[S(n, m) = sum_{i = 1}^m varphi(ni) ]

考虑如何将 (varphi(ni)) 给拆开，因为 (varphi(x)) 只有每个质因子 (p) 第一次出现的时候才会是 (p - 1) ，其他时候都直接乘 (p) 。

我们设 $n = prod_i p_i^{a_i} $ ，那么令 (p = prod_i p_i^{a_i - 1}, q = prod_i p_i) 那么有 (pq = n) 。

有一些预备知识可以看看 这篇博客 ，讲了一下如何证明后面的一些结论。

假设你前置的欧拉函数知识都会了，那么我接下来就把 yyb 的博客搬过来啦。

[egin{aligned} S(n,m) &= p sum_{i = 1}^m varphi(qi)\ &= psum_{i = 1}^m varphi(q) varphi(frac{i}{gcd(i,q)}) gcd(i,q)\ &= p sum_{i = 1}^m varphi(frac q{gcd(i,q)}) varphi(i) gcd(i,q)\ &= p sum_{i = 1}^m varphi(frac q{gcd(i,q)}) varphi(i) sum_{d|gcd(i,q)} varphi(d)\ &= p sum_{i = 1}^m varphi(i) sum_{d|i,d|q}varphi(frac{q}{d})\ &= p sum_{d | q} varphi(frac qd) sum_{i=1}^{frac md} varphi(id)\ &= p sum_{d | q} varphi(frac qd) S(d, lfloor frac{m}{d} floor)\ end{aligned} ]

其中有几步比较 ( ext{tricky}) ，就比如 (d | gcd(i, q)) 变回 (d | i, d | q) 然后枚举 (d | q) ，算答案。

最后我们利用这个式子直接记忆化处理就行了，对于 (n = 1) 的时候用杜教筛处理边界就行了。

复杂度有个比较松的上界 (mathcal O(n sqrt m + m^frac 2 3)) ，但好像太松了。。我极限只需要 (0.6s) 。（还是用的 std :: map 。。）

代码

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
#define pb push_back

using namespace std;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

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

void File() {
#ifdef zjp_shadow
	freopen ("3512.in", "r", stdin);
	freopen ("3512.out", "w", stdout);
#endif
}

const int Lim = 1e5, N = Lim + 1e3, Mod = 1e9 + 7;

int prime[N], pcnt, phi[N], sumphi[N], minp[N];
bitset<N> is_prime;

void Linear_Sieve(int maxn) {
	is_prime.set(); is_prime[0] = is_prime[1] = false;
	phi[1] = sumphi[1] = 1;
	For (i, 2, maxn) {
		if (is_prime[i]) prime[++ pcnt] = i, phi[i] = i - 1, minp[i] = i;
		for (int j = 1, res; (res = i * prime[j]) <= maxn && j <= pcnt; ++ j) {
			is_prime[res] = false; minp[res] = prime[j];
			if (i % prime[j]) phi[res] = phi[i] * (prime[j] - 1);
			else { phi[res] = phi[i] * prime[j]; break; }
		}
		sumphi[i] = (sumphi[i - 1] + phi[i]) % Mod;
	}
}

map<int, int> M, val[N];

#define Out(a, b) if (a) return b;

int Phi(int n) {
	Out(n <= Lim, sumphi[n]); Out(M[n], M[n]);
	int res = 1ll * n * (n + 1) / 2 % Mod;
	for (int i = 2, ni; i <= n; i = ni + 1)
		res = (res + Mod - ((ni = n / (n / i)) - i + 1ll) * Phi(n / i)) % Mod;
	return M[n] = res;
}

int S(int n, int m) {
	Out(!m, 0); Out(n == 1, Phi(m)); Out(m == 1, phi[n]); Out(val[n][m], val[n][m]);
	int res = 0, p = 1, q = 1, tmp = n; vector<int> fac;
	while (tmp > 1) {
		int x = minp[tmp]; q *= x; tmp /= x; fac.pb(x);
		while (!(tmp % x)) p *= x, tmp /= x;
	}
	Rep (i, 1 << int(fac.size())) {
		int d = 1;
		Rep (j, fac.size()) if (i >> j & 1) d *= fac[j];
		res = (res + 1ll * phi[q / d] * S(d, m / d)) % Mod;
	}
	return val[n][m] = 1ll * res * p % Mod;
}

int main () {

	File();

	Linear_Sieve(Lim);

	int n = read(), m = read(), ans = 0;
	For (i, 1, n) 
		ans = (ans + S(i, m)) % Mod;
	printf ("%d
", ans);

	return 0;

}

BZOJ 3560: DZY Loves Math V

题意

给你 (n) 个正整数 (a_1, a_2, cdots, a_n) ，求

[sum_{i_1 | a_1} sum_{i_2 | a_2} cdots sum_{i_n | a_n} varphi(i_1 i_2 cdots i_n) pmod {10^9 + 7} ]

(1 le n le 10^5，1 le a_i le 10^7)

题解

终于不是那么神的题啊。。。

又发现答案又是一个积性函数。

然后显然对于每个质因子 (p) 单独考虑，假设其在 (a_i) 中出现 (b_i) 次。

那么答案其实就是

[ans = prod_p sum_{i_1 = 0}^{b_1} sum_{i_2 = 0}^{b_2} dots sum_{i_n = 0}^{b_n} varphi(p^{sum_{j = 1}^n i_j}) ]

考虑当 (k > 0) 时有 (varphi(p^k) = p^k imes displaystyle frac {p - 1} p) ，对于 (k = 0) 时则为 (varphi(1) = 1) 。

那么拆一下就变成

[ans = prod_p (((sum_{i_1 = 0}^{b_1} sum_{i_2 = 0}^{b_2} dots sum_{i_n = 0}^{b_n} p^{sum_{j = 1}^n i_j}) - 1) imes frac{p - 1}p + 1) ]

然后对于直接拆回去就行了。

[ans = prod_p ((prod_{i = 1}^n sum_{j = 0}^{b_i} p^j) - 1) imes frac{p - 1}p + 1) ]

然后分解质因数逐个乘起来就行了，复杂度暴力点是 (mathcal O(n sqrt {a_i})) 的。

代码

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
#define mp make_pair
#define fir first
#define sec second

using namespace std;

typedef pair<int, int> PII;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

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

void File() {
#ifdef zjp_shadow
	freopen ("3560.in", "r", stdin);
	freopen ("3560.out", "w", stdout);
#endif
}

const int Lim = 1e7, N = Lim + 1e3, Mod = 1e9 + 7;

inline int fpm(int x, int power) {
	int res = 1;
	for (; power; power >>= 1, x = 1ll * x * x % Mod)
		if (power & 1) res = 1ll * res * x % Mod;
	return res;
}

PII fac[N]; int ftot, inv[N];

void ReSolve(int x) {
	for (int i = 2; i * i <= x; ++ i) if (!(x % i)) {
		int power = 0;
		while (!(x % i)) x /= i, ++ power;
		fac[++ ftot] = mp(i, power);
	}
	if (x > 1) fac[++ ftot] = mp(x, 1);
}

int sump[210];

int Solve(int l, int r, int p) {
	int res = 1, cur = 0, phip = 1; sump[0] = 1;
	For (i, l, r) {
		for (; cur < fac[i].sec; ++ cur) {
			phip = 1ll * phip * p % Mod;
			sump[cur + 1] = (sump[cur] + phip) % Mod;
		}
		res = 1ll * res * sump[fac[i].sec] % Mod;
	}
	return 1ll * (res - 1) * (p - 1) % Mod * fpm(p, Mod - 2) % Mod + 1;
}

int main () {

	File();

	int n = read();
	For (i, 1, n) ReSolve(read());
	sort(fac + 1, fac + ftot + 1);

	int ans = 1;
	for (int i = 1, j; i <= ftot; i = j + 1) {
		for (j = i; j < ftot && fac[j + 1].fir == fac[i].fir; ++ j);
		ans = 1ll * ans * Solve(i, j, fac[i].fir) % Mod;
	}
	printf ("%d
", ans);

	return 0;

}

BZOJ 3561: DZY Loves Math VI

题意

给定整数 (n, m) 求

[sum_{i = 1}^{n} sum_{j = 1}^{m} mathrm{lcm}(i, j)^{gcd(i, j)} pmod {10^9 + 7} ]

(n, m le 5 imes 10^5)

题解

傻吊题。

[ans = sum_{d = 1}^{n} d^d sum_{x = 1}^{min(lfloor frac n d floor, lfloor frac m d floor)} mu(x) x^{2d} (sum_{i = 1}^{lfloor frac n d floor} i^d) (sum_{j = 1}^{lfloor frac m d floor} j^d) pmod {10^9 + 7} ]

考虑从小到大枚举 (d) ，然后每次把 (displaystyle 1 sim lfloor frac n d floor) 的 (i^{d - 1}) 乘上 (i) 更新到 (i^d) 就行了，然后顺便算下前缀和就行了。

复杂度是调和级数 (mathcal O(n log n)) 的。

BZOJ 3568: DZY Loves Math VII

题意

已知 (mu(n)) ，求第 (k) 小的 (n) 。( ext{I.} |mu(n)| le 1,k le 10^8)

已知 (varphi(n)) ，求第 (k) 小的 (n) 。( ext{II.} varphi(n) le 10^{10}, kle 1000)

已知 (d(n)) ，求第 (k) 小的 (n) 。( ext{III.} d(N) le 10^7,K le 50)

题解

(mu(n))

第一问考虑二分 (n) ，对于 (sum_{i = 1}^{n} mu(i)) 杜教筛即可。

然后这只能算出和，不能确定 ({-1, 0, 1}) 分别有多少，我们再算一下 (sum_{i = 1}^{n} |mu(i)|) 的大小就行了，就能先确定 (0) 的个数，然后列方程确定 ({-1, 1}) 的个数了。

如何算呢？容斥一下就行了。

[sum_{i = 1}^{n} |mu(i)| = n + sum_{i = 2}^{lfloor sqrt n floor} mu(i) lfloor frac n {i^2} floor ]

大概意思就是考虑每个含有平方因子的 (a) ，都不能被算进去，我们要减掉。

那么利用 (mu(x)) 的容斥性质，在一个质因子的时候算 (-1) 减掉，两个加上 (cdots)

(varphi(n))

我们利用 (varphi(n)) 的计算式。

[varphi(n) = prod_p p^k frac{p - 1}{p} ]

大力分解 (varphi(n)) ，然后直接算所有合法解，小于 (10^6) 暴力枚举试除，大于 (10^6) 用 Miller_Rabin 判断是否合法就行了，因为 (> 10^6) 的质因子不可能超过一个。

看起来很暴力，其实可以通过。。

(d(n))

反质数加强版。。。

利用 (d(n)) 的计算式

[d(n) = prod_p (k + 1) ]

大力分解一波，(K = 1) 时候指数应该是单调不增的，不然一定可以搞到一个更优解，然后大小判断用 (ln) 。

那么我们基于这个策略进行调整，然后剪剪枝，还需要一个高精度。

这题难点应该就在这问，实现起来十分麻烦。而且剪枝看起来十分不靠谱。。。