BZOJ4407 于神之怒加强版

题解

非常裸的莫比乌斯反演。

但是反演完还需要快速计算一个积性函数（我直接用$nlogn$卷积被TLE了

推荐一个博客

我也不想再写一遍了

代码

#include<cstring>
#include<cstdio>
#include<algorithm>
#define ll long long
#define rd read()
using namespace std;
 
const ll mod = 1e9 + 7;
 
const int N = 5e6;
 
int T, n, m, k;
int mu[N], pri[N], tot, vis[N], low[N];
ll sum[N];
 
int read() {
    int X = 0, p = 1; char c = getchar();
    for(; c > '9' || c < '0'; c = getchar()) if(c == '-') p = -1;
    for(; c >= '0' && c <= '9'; c = getchar()) X = X * 10 + c - '0';
    return X * p;
}
 
ll fpow(ll a, ll b) {
    ll re = 1;
    for(; b; b >>= 1, a = a * a % mod) if(b & 1) re = re * a % mod;
    return re;
}
 
void init() {
    mu[1] = 1;
    sum[1] = 1;
    for(int i = 2; i < N; ++i) {
        if(!vis[i]) mu[i] = -1, pri[++tot] = i, sum[i] = fpow(i, k) - 1, low[i] = i;
        for(int j = 1; j <= tot && i * pri[j] < N; ++j) {
            vis[i * pri[j]] = 1;
            if(i % pri[j] == 0) {
                low[i * pri[j]] = low[i] * pri[j];
                mu[i * pri[j]] = 0;
                if(low[i] == i) sum[i * pri[j]] = sum[i] * fpow(pri[j], k) % mod;
                else sum[i * pri[j]] =  sum[i / low[i]] * sum[pri[j] * low[i]] % mod;
                break;
            }
            low[i * pri[j]] = pri[j];
            mu[i * pri[j]] = -mu[i];
            sum[i * pri[j]] = sum[i] * sum[pri[j]] % mod;
        }
    }
    //for(ll i = 1; i < N; ++i) {
    //  ll tmp = fpow(i, k);
    //  for(ll j = 1; j < N && i * j < N; ++j) sum[i * j] = (sum[i * j] + tmp * mu[j]) % mod;
    //}
    for(int i = 1; i < N; ++i) sum[i] = (sum[i - 1] + sum[i]) % mod;
}
 
ll cal() {
    ll ans = 0;
    for(int i = 1, j = 1; i <= min(n, m); i = j + 1) {
        j = n / (n / i);
        j = min(j, m / (m / i));
        ans += (sum[j] - sum[i - 1]) % mod * (n / i) % mod * (m / i) % mod;
        ans %= mod;
    }
    ans = (ans % mod + mod) % mod;
    return ans;
}
 
int main()
{
    T = rd; k = rd;
    init();
    for(; T; T--) {
        n = rd; m = rd;
        printf("%lld
", cal());
    }
}