[Cqoi 2017] bzoj4815 小Q的表格 [莫比乌斯反演]

Description:
自己看题。
Solution:
这个$f$的形式和$gcd$很像，事实上每次修改$(a,b)$会影响到$gcd(i,j)==gcd(a,b)$。
所以我们只用保留$gcd(i,i)$即可。每次修改$(a,b)$等价于修改$gcd(a,b)$。
然后维护$S(d)$，表示$f(d,d)$的前缀和。
那么得出答案即是
$\sum_{i=1}^{k}{f(i)S([\frac{k}{i}])}$
$S(n)=\sum_{i=1}^{n}{\phi(i)*i^2}$

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int maxn = 4e6 + 5, P = 1e9 + 7;
int n, m;
int mark[maxn], p[maxn], l[maxn], r[maxn], b[maxn];
ll phi[maxn], s[maxn], S[maxn], add[maxn], a[maxn];
int gcd(int a, int b) {
    return !b ? a : gcd(b, a % b);
}
ll power(ll x, ll t) {
    ll ret = 1;
    for(; t; t >>= 1, x = x * x % P) {
        if(t & 1) {
            ret = ret * x % P;
        }
    }
    return ret;
}
void ini() {
    phi[1] = 1;
    for(int i = 2; i < maxn; ++i) {
        if(!mark[i]) {
            p[++p[0]] = i;
            phi[i] = i - 1;
        }
        for(int j = 1; j <= p[0] && i * p[j] < maxn; ++j) {
            mark[i * p[j]] = 1;
            if(i % p[j] == 0) {
                phi[i * p[j]] = phi[i] * p[j];
                break;
            }
            phi[i * p[j]] = phi[i] * phi[p[j]];
        }
    }
    for(int i = 1; i < maxn; ++i) {
        s[i] = (s[i - 1] + phi[i] * i % P * i % P) % P;
    }
}
int main() {
    ini();
    scanf("%d%d", &m, &n);
    int B = sqrt(n);
    for(int i = 1; i <= n; ++i) {
        S[i] = (S[i - 1] + (long long) i * i % P) % P;
        a[i] = (long long) i * i % P;
        b[i] = (i - 1) / B + 1;
        if(!l[b[i]]) {
            l[b[i]] = i;
        }
        r[b[i]] = i;
    }
    while(m--) {
        int A, B, k;
        ll x, ans = 0;
        scanf("%d%d%lld%d", &A, &B, &x, &k);
        int d = gcd(A, B);
        x %= P;
        x = x * d % P * d % P * power((long long) A * B % P, P - 2) % P;
        for(int i = d; i <= r[b[d]]; ++i) {
            S[i] = (S[i] + x - a[d] + P) % P;
        }
        for(int i = b[d] + 1; i <= b[n]; ++i) {
            add[i] = (add[i] + x - a[d] + P) % P;
        }
        a[d] = x;
        for(int i = 1, j; i <= k; i = j + 1) {
            j = k / (k / i);
            ans = (ans + (S[j] + add[b[j]] - S[i - 1] - add[b[i - 1]] + 2 * P) % P * s[k / i] % P) % P;
        }
        printf("%lld\n", ans);
    }       
    return 0;
}