Luogu 4449 于神之怒加强版

挺套路的题，然而一开始还是想错了……

$\sum_{i = 1}^{n}\sum_{j = 1}^{m}gcd(i, j) ^ {k} = \sum_{T = 1}^{min(n, m)}\left \lfloor \frac{n}{T} \right \rfloor \left \lfloor \frac{m}{T} \right \rfloor\sum_{d | T}\mu (\frac{T}{d}) * d^{k}$

我们设$h(i) =\sum_{d | T}\mu (\frac{T}{d}) * d^{k} $， 那么原式化为$ \sum_{T = 1}^{min(n, m)}\left \lfloor \frac{n}{T} \right \rfloor \left \lfloor \frac{m}{T} \right \rfloor h(i)$

我们做出$h(i)$的前缀和之后就可以整除分块了。

发现$h(i)$是一个积性函数，可以线性筛，具体筛法如下：

1、$h(1) = 1$

2、$h(p) = p^{k} - 1$($p$是一个质数)

3、考虑线性筛中的过程，$i * pri_{j}$被它的最小因子也就是$pri_{j}$筛掉，那么当$i$不是$pri_{j}$的倍数的时候，也就是说$h(i * pri_{j}) = h(i) * h(pri_{j})$。

　　而当$pri_{j} | i$的时候，考虑分两步：（$ i = \prod_{j = 1}^{m}p_{j}^{c_{j}}$）

　　记$low_{i} =$ $i$的最小质因子被$i$包含的最大幂次积（即$p_{1} ^ {c_{1}}$）

　　a、我们假设$low_{i}$不等于$i$，这句话等价于$i$不是一个质数的倍数，所以$\frac{i}{low_{i}}$ 与 $pri_{j} * low_{i}$互质，那么直接乘上就好了。

　　b、当$low_{i}$等于$i$的时候，我们考虑一下$h(i)$函数的定义，展开式子之后发现$h(i) = pri_{j} ^ {mk} - pri_{j} ^ {(m - 1)k}$，而$h(i  * pri_{j}) = pri_{j} ^ {(m + 1)k} - pri_{j} ^ {mk}$，那么就相当于$h(i * pri_{j}) = h(i) * pri_{j}^{k}$。

时间复杂度$O(maxNlogmaxN + T\sqrt{n})$。

大时限题还感觉跑得挺快的。

Code：

#include <cstdio>
#include <cstring>
using namespace std;
typedef long long ll;

const int N = 5e6 + 5;
const ll P = 1e9 + 7;

int testCase, pCnt = 0, pri[N];
ll k, h[N], low[N], sum[N];
bool np[N];

template <typename T>
inline void read(T &X) {
    X = 0;
    char ch = 0;
    T op = 1;
    for(; ch > '9' || ch < '0'; ch = getchar())
        if(ch == '-') op = -1;
    for(; ch >= '0' && ch <= '9'; ch = getchar())
        X = (X << 3) + (X << 1) + ch - 48;
    X *= op;
}

inline ll pow(ll a, ll b) {
    ll res = 1LL;
    for(; b > 0; b >>= 1) {
        if(b & 1) res = res * a % P;
        a = a * a % P;
    }
    return res;
}

inline void sieve() {
    h[1] = 1LL;
    for(int i = 2; i < N; i++) {
        if(!np[i]) {
            pri[++pCnt] = i;
            low[i] = (ll)i;
            h[i] =
             (pow(i, k) - 1LL + P) % P;
        }
        for(int j = 1; j <= pCnt && pri[j] * i < N; j++) {
            np[i * pri[j]] = 1;
            if(i % pri[j] == 0) {
                low[i * pri[j]] = low[i] * pri[j];
                if(low[i] == i) h[i * pri[j]] = h[i] * pow(pri[j], k) % P;
                else h[i * pri[j]] = h[i / low[i]] * h[low[i] * pri[j]] % P;
                break;
            }
            low[i * pri[j]] = pri[j];
            h[i * pri[j]] = h[i] * h[pri[j]] % P;
        }
    }
    
    
/*    for(int i = 1; i < 20; i++)
        printf("%lld ", phi[i]);
    printf("\n");   
    for(int i = 1; i < 30; i++)
        printf("%lld ", h[i]);
    printf("\n");   */
    
/*    for(int i = 1; i <= 20; i++)
        printf("%lld ", h[i] % P);   */
    
    for(int i = 1; i < N; i++) 
        sum[i] = (sum[i - 1] + h[i] % P + P) % P;
}

inline ll min(ll x, ll y) {
    return x > y ? y : x;
}

int main() {
    read(testCase), read(k);
    sieve();
    for(ll n, m; testCase--; ) {
        read(n), read(m);
        ll ans = 0LL, rep = min(n, m);
        for(ll l = 1, r; l <= rep; l = r + 1) {
            r = min((n / (n / l)), (m / (m / l)));
            ans = (ans + (sum[r] - sum[l - 1] + P) % P * (n / l) % P * (m / l) % P) % P;
        }
        printf("%lld\n", ans);
    } 
    return 0;
}