# [BZOJ4407] 于神之怒加强版

Portal

$Ans = \sum_{i = 1}^{n} \sum_{j = 1}^{m} (i, j) ^{k} \pmod {1e9 + 7}\\$

$Ans = \sum_{i}\sum_{j}\sum_{l = 1} ^ {min(n, m)} l ^ {k} ~\cdot~ [(i, j) == l] \\ = \sum_{l = 1} ^ {min(n, m)} l ^ {k} \sum_{i}\sum_{j} [(\frac{i}{l}, \frac{j}{l})= 1] \\ = \sum_{l = 1} ^ {min(n, m)} l ^ {k} \sum_{i}^{\frac{n}{l}}\sum_{j}^{\frac{m}{l}} [(i,j)= 1] \\ = \sum_{l = 1} ^ {min(n, m)} l ^ {k} \sum_{i}\sum_{j} \sum_{d} [d | i][d | j] \mu(d) \\ = \sum_{l}l^k \sum_{d} \mu(d)\frac{n}{ld}\frac{m}{ld} \\$

$Ans = \sum_{l}\sum_{d} l^k\mu(d)\frac{n}{ld}\frac{m}{ld}\\ =\sum_{c}\sum_{p | c} p ^ k \frac{}{} \mu(\frac{c}{p}) \frac{}{} \frac{n}{c} \frac{m}{c}$

$(1)$ 如果 i % prime[j] != 0, F[i * prime[j]] = F[i] * F[prime[j]]

### Code

#include<bits/stdc++.h>
using namespace std;
#define rep(i, a, b) for(int i = (a), i##_end_ = (b); i <= i##_end_; ++i)
#define drep(i, a, b) for(int i = (a), i##_end_ = (b); i >= i##_end_; --i)
#define clar(a, b) memset((a), (b), sizeof(a))
#define debug(...) fprintf(stderr, __VA_ARGS__)
typedef long long LL;
typedef long double LD;
char ch = getchar();
int x = 0, flag = 1;
for (;!isdigit(ch); ch = getchar()) if (ch == '-') flag *= -1;
for (;isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
return x * flag;
}
void write(int x) {
if (x < 0) putchar('-'), x = -x;
if (x >= 10) write(x / 10);
putchar(x % 10 + 48);
}

const int Mod = 1e9 + 7;
const int Maxn = 5e6 + 9;
int fpm(int base, int tims) {
int r = 1; base %= Mod;
while (tims) {
if (tims & 1) r = 1ll * base * r % Mod;
base = 1ll * base * base % Mod;
tims >>= 1;
}
return r;
}

int T, k;

int prime[Maxn], isnprime[Maxn], mu[Maxn], tot, F[Maxn], low[Maxn], pk[Maxn];
int sumfixF[Maxn];
inline int getPk(int u) { return pk[u] ? pk[u] : (pk[u] = fpm(u, k)); }
void linearSieve() {
mu[1] = 1; F[1] = 1;
rep (i, 2, Maxn - 1) {
if (!isnprime[i]) prime[++tot] = i, mu[i] = -1, F[i] = (getPk(i) - 1ll + Mod) % Mod, low[i] = i;
for (int k, j = 1; j <= tot && (k = prime[j] * i) < Maxn; ++j) {
isnprime[k] = 1;
if (i % prime[j] == 0) {
low[k] = low[i] * prime[j], mu[k] = 0;
if (low[k] == k) F[k] = 1ll * F[i] * getPk(prime[j]) % Mod;
/**/			else F[k] = 1ll * F[k / low[k]] * F[low[k]] % Mod;
break;
} else {
mu[k] = mu[i] * mu[prime[j]];
F[k] = 1ll * F[i] * F[prime[j]] % Mod;
low[k] = prime[j];
}
}
}

rep (i, 1, Maxn - 1) sumfixF[i] = (1ll * sumfixF[i - 1] + F[i]) % Mod;

}

vector <pair<int, int> > s;
void init() {
rep (i, 1, T) {
s.push_back(make_pair(u, v));
}

linearSieve();
}

void solve() {
rep (i, 1, T) {
int n = s[i - 1].first, m = s[i - 1].second;
int Limit = min(n, m); LL ans = 0;
for (int r, l = 1; l <= Limit; l = r + 1) {
r = min(min(n / (n / l), m / (m / l)), Limit);
(ans += 1ll * (1ll * sumfixF[r] - sumfixF[l - 1] + Mod) % Mod * (n / l) % Mod * (m / l) % Mod) %= Mod;
}
printf("%d\n", ans);
}
}

int main() {
freopen("BZOJ4407.in", "r", stdin);
freopen("BZOJ4407.out", "w", stdout);

init();
solve();

#ifdef Qrsikno
debug("\nRunning time: %.3lf(s)\n", clock() * 1.0 / CLOCKS_PER_SEC);
#endif
return 0;
}


posted @ 2018-12-29 17:33  Qrsikno  阅读(110)  评论(0编辑  收藏  举报