YY 的 GCD
题目大意
求下列式子的值:
∑ i = 1 n ∑ j = 1 m [ gcd { i , j } ∈ P ] ( P 为素数集合 ) \sum_{i = 1}^{n} \sum_{j = 1}^{m} [\gcd\{i, j\} \in P]\ (P\ \text{为素数集合}) i=1∑nj=1∑m[gcd{i,j}∈P] (P 为素数集合)
1 0 4 10^4 104 组询问, n , m ≤ 1 0 7 n, m \le 10^7 n,m≤107。
思路分析
将原式变形:
∑ p ∈ P ∑ i = 1 ⌊ n p ⌋ ∑ j = 1 ⌊ m p ⌋ [ gcd { i , j } = 1 ] \sum_{p \in P} \sum_{i = 1}^{\lfloor \frac{n}{p} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{p} \rfloor} [\gcd\{i, j\} = 1] p∈P∑i=1∑⌊pn⌋j=1∑⌊pm⌋[gcd{i,j}=1]
∑ p ∈ P ∑ i = 1 ⌊ n p ⌋ ∑ j = 1 ⌊ m p ⌋ ∑ d ∣ i , d ∣ j μ ( d ) \sum_{p \in P} \sum_{i = 1}^{\lfloor \frac{n}{p} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{p} \rfloor} \sum_{d | i, d | j} \mu(d) p∈P∑i=1∑⌊pn⌋j=1∑⌊pm⌋d∣i,d∣j∑μ(d)
∑ p ∈ P ∑ d = 1 min { ⌊ n p ⌋ , ⌊ n p ⌋ } μ ( d ) ∑ i = 1 ⌊ n p d ⌋ ∑ j = 1 ⌊ m p d ⌋ 1 \sum_{p \in P} \sum_{d = 1}^{\min\{\lfloor \frac{n}{p} \rfloor, \lfloor \frac{n}{p} \rfloor\}} \mu(d) \sum_{i = 1}^{\lfloor \frac{n}{pd} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{pd} \rfloor} 1 p∈P∑d=1∑min{⌊pn⌋,⌊pn⌋}μ(d)i=1∑⌊pdn⌋j=1∑⌊pdm⌋1
∑ p ∈ P ∑ d = 1 min { ⌊ n p ⌋ , ⌊ n p ⌋ } μ ( d ) ⌊ n p d ⌋ ⌊ m p d ⌋ \sum_{p \in P} \sum_{d = 1}^{\min\{\lfloor \frac{n}{p} \rfloor, \lfloor \frac{n}{p} \rfloor\}} \mu(d) \lfloor \frac{n}{pd} \rfloor \lfloor \frac{m}{pd} \rfloor p∈P∑d=1∑min{⌊pn⌋,⌊pn⌋}μ(d)⌊pdn⌋⌊pdm⌋
∑ T = 1 min { n , m } ∑ p ∣ T , p ∈ P μ ( T p ) ⌊ n T ⌋ ⌊ m T ⌋ \sum_{T = 1}^{\min\{n, m\}} \sum_{p | T, p \in P} \mu(\frac{T}{p}) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor T=1∑min{n,m}p∣T,p∈P∑μ(pT)⌊Tn⌋⌊Tm⌋
∑ T = 1 min { n , m } ⌊ n T ⌋ ⌊ m T ⌋ ∑ p ∣ T , p ∈ P μ ( T p ) \sum_{T = 1}^{\min\{n, m\}} \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \color{red}{\sum_{p | T, p \in P} \mu(\frac{T} {p})} T=1∑min{n,m}⌊Tn⌋⌊Tm⌋p∣T,p∈P∑μ(pT)
红色部分可以使用线性筛算出 μ \mu μ,然后枚举每个素数并更新它的倍数算出,时间复杂度 O ( n log log n ) O(n \log \log n) O(nloglogn)。黑色部分可以使用数论分块的技巧求出。
代码实现
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long llong;
const int maxn = 1e7;
int T, n, m, p[maxn + 3];
int mrk[maxn + 3], miu[maxn + 3], sum[maxn + 3];
int main() {
miu[1] = 1;
for (int i = 2; i <= maxn; i++) {
if (!mrk[i]) p[++m] = i, miu[i] = -1;
for (int j = 1, t; j <= m && (t = i * p[j]) <= maxn; j++) {
mrk[t] = 1;
if (i % p[j] == 0) break;
miu[t] = -miu[i];
}
}
for (int i = 1; i <= m; i++) {
for (int j = 1, t; (t = p[i] * j) <= maxn; j++) {
sum[t] += miu[j];
}
}
for (int i = 1; i <= maxn; i++) {
sum[i] += sum[i - 1];
}
scanf("%d", &T);
while (T--) {
scanf("%d %d", &n, &m);
if (n > m) swap(n, m);
llong ans = 0;
for (int i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i), m / (m / i));
ans += 1ll * (n / i) * (m / i) * (sum[j] - sum[i - 1]);
}
printf("%lld\n", ans);
}
return 0;
}
NOI 2016 循环之美
题目大意
求下列式子的值:
∑ i = 1 n ∑ j = 1 m [ gcd { i , j } = 1 ] [ gcd { j , k } = 1 ] \sum_{i = 1}^{n} \sum_{j = 1}^{m} [\gcd\{i, j\} = 1] [\gcd\{j, k\} = 1] i=1∑nj=1∑m[gcd{i,j}=1][gcd{j,k}=1]
n , m ≤ 1 0 9 , k ≤ 2000 n, m \le 10^9, k \le 2000 n,m≤109,k≤2000。
思路分析
将原式变形:
∑ i = 1 m [ gcd { i , k } = 1 ] ∑ j = 1 n [ gcd { i , j } = 1 ] \sum_{i = 1}^{m} [\gcd\{i, k\} = 1] \sum_{j = 1}^{n} [\gcd\{i, j\} = 1] i=1∑m[gcd{i,k}=1]j=1∑n[gcd{i,j}=1]
∑ i = 1 m [ gcd { i , k } = 1 ] ∑ j = 1 n ∑ d ∣ i , j μ ( d ) \sum_{i = 1}^{m} [\gcd\{i, k\} = 1] \sum_{j = 1}^{n} \sum_{d | i, j} \mu(d) i=1∑m[gcd{i,k}=1]j=1∑nd∣i,j∑μ(d)
∑ d = 1 min { n , m } μ ( d ) ∑ i = 1 ⌊ m d ⌋ [ gcd { i ⋅ d , k } = 1 ] ∑ j = 1 ⌊ n d ⌋ 1 \sum_{d = 1}^{\min\{n, m\}} \mu(d) \sum_{i = 1}^{\lfloor \frac{m}{d} \rfloor}[\gcd\{i \cdot d, k\} = 1] \sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor} 1 d=1∑min{n,m}μ(d)i=1∑⌊dm⌋[gcd{i⋅d,k}=1]j=1∑⌊dn⌋1
∑ d = 1 min { n , m } μ ( d ) ⋅ [ gcd { d , k } = 1 ] ⋅ ⌊ n d ⌋ ∑ i = 1 ⌊ m d ⌋ [ gcd { i , k } = 1 ] \sum_{d = 1}^{\min\{n, m\}} \color{blue}{\mu(d) \cdot [\gcd\{d, k\} = 1]} \cdot \lfloor \frac{n}{d} \rfloor\sum_{i = 1}^{\lfloor \frac{m}{d} \rfloor}\color{red}{[\gcd\{i, k\} = 1]} d=1∑min{n,m}μ(d)⋅[gcd{d,k}=1]⋅⌊dn⌋i=1∑⌊dm⌋[gcd{i,k}=1]
现在瓶颈就变成了计算红色部分和蓝色部分的前缀和。
注意到 k ≤ 2000 k \le 2000 k≤2000,于是红色部分直接可以递推计算 ≤ k \le k ≤k 的部分,然后在乘上计数的次数即可。
对于蓝色部分,令 S ( n , k ) = ∑ i = 1 n [ gcd ( i , k ) = 1 ] ⋅ μ ( i ) S(n, k) = \sum_{i = 1}^{n} [\gcd(i, k) = 1] \cdot \mu(i) S(n,k)=∑i=1n[gcd(i,k)=1]⋅μ(i)。对原式进行变形,得到原式 = ∑ d ∣ k μ ( d ) ∑ i = 1 ⌊ n d ⌋ μ ( i ⋅ d ) = \sum_{d | k} \mu(d) \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(i \cdot d) =∑d∣kμ(d)∑i=1⌊dn⌋μ(i⋅d)。发现如果 gcd ( i , d ) ≠ 1 \gcd(i, d) \neq 1 gcd(i,d)̸=1,那么 μ ( i ⋅ d ) = 0 \mu(i \cdot d) = 0 μ(i⋅d)=0,对答案没有影响。于是可以继续化简,得到 S ( n , k ) = ∑ d ∣ k μ ( d ) 2 S ( ⌊ n d ⌋ , d ) S(n, k) = \sum_{d | k} \mu(d) ^ 2 S(\lfloor \frac{n}{d} \rfloor, d) S(n,k)=∑d∣kμ(d)2S(⌊dn⌋,d)。直接记忆化搜索 + + + 杜教筛即可。
代码实现
#include <cstdio>
#include <map>
#include <algorithm>
using namespace std;
typedef pair<int, int> pii;
const int maxm = 1e7, maxk = 2000;
bool mrk[maxm + 3];
int n, m, k, cnt, p[maxm + 3], miu[maxm + 3], sum[maxm + 3], f[maxk + 3];
map<pii, int> mp;
int gcd(int a, int b) {
return b ? gcd(b, a % b) : a;
}
int F(int n) {
return f[n % k] + n / k * f[k];
}
int S(int n, int k) {
if ((k == 1 && n <= maxm) || !n) {
return sum[n];
}
pii x = pii(n, k);
if (mp.find(x) != mp.end()) {
return mp[x];
}
int res;
if (k == 1) {
res = 1;
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
res -= (r - l + 1) * S(n / l, 1);
}
} else {
res = 0;
for (int i = 1, j; i * i <= k; i++) {
if (k % i == 0) {
j = k / i;
if (miu[i]) res += S(n / i, i);
if (miu[j] && i * i != k) res += S(n / j, j);
}
}
}
return mp[x] = res;
}
int main() {
#ifdef LOCAL
freopen("input.md", "r", stdin);
freopen("output.md", "w", stdout);
#endif
miu[1] = 1;
for (int i = 2; i <= maxm; i++) {
if (!mrk[i]) {
p[++cnt] = i;
miu[i] = -1;
}
for (int j = 1, t; j <= cnt && (t = i * p[j]) <= maxm; j++) {
mrk[t] = 1;
if (i % p[j] == 0) break;
miu[t] = -miu[i];
}
}
for (int i = 1; i <= maxm; i++) {
sum[i] = sum[i - 1] + miu[i];
}
scanf("%d %d %d", &n, &m, &k);
for (int i = 1; i <= k; i++) {
f[i] = f[i - 1] + (gcd(i, k) == 1);
}
int x = min(n, m);
long long ans = 0;
for (int l = 1, r; l <= x; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans += 1ll * (n / l) * F(m / l) * (S(r, k) - S(l - 1, k));
}
printf("%lld\n", ans);
return 0;
}
UOJ 62 怎样跑的更快
题目大意
神仙题,参考官方题解。
思路分析
代码实现
#include <cstdio>
const int maxn = 1e5, mod = 998244353;
int n, c, d, q, f_r[maxn + 3], f_z[maxn + 3], hx[maxn + 3];
int qpow(int a, int b) {
int c = 1;
for (; b; b >>= 1, a = 1ll * a * a % mod) {
if (b & 1) c = 1ll * a * c % mod;
}
return c;
}
void add(int &a, int b) {
a += b;
a < mod ? NULL : a -= mod;
}
void sub(int &a, int b) {
a -= b;
a < 0 ? a += mod : NULL;
}
void solve(int n, int f[], int type) {
if (!type) {
for (int i = 1; i <= n; i++) {
for (int j = 2 * i; j <= n; j += i) {
sub(f[j], f[i]);
}
}
} else {
for (int i = n; i; i--) {
for (int j = 2 * i; j <= n; j += i) {
sub(f[i], f[j]);
}
}
}
}
int main() {
#ifdef LOCAL
freopen("input.md", "r", stdin);
freopen("output.md", "w", stdout);
#endif
scanf("%d %d %d %d", &n, &c, &d, &q);
c %= mod - 1, d %= mod - 1;
c = (c - d + mod - 1) % (mod - 1);
for (int i = 1; i <= n; i++) {
f_r[i] = qpow(i, c);
}
solve(n, f_r, 0);
while (q--) {
for (int i = 1; i <= n; i++) {
scanf("%d", &f_z[i]);
f_z[i] = 1ll * f_z[i] * qpow(i, mod - 1 - d) % mod;
}
solve(n, f_z, 0);
bool flag = true;
for (int i = 1; i <= n; i++) {
if (f_r[i] > 0) {
hx[i] = 1ll * f_z[i] * qpow(f_r[i], mod - 2) % mod;
} else if (f_z[i] > 0) {
flag = false;
break;
}
}
if (!flag) {
puts("-1");
continue;
}
solve(n, hx, 1);
for (int i = 1; i <= n; i++) {
printf("%d%c", 1ll * hx[i] * qpow(i, mod - 1 - d) % mod, " \n"[i == n]);
}
}
return 0;
}