CQOI 2017 小Q的表格 - 不一样的暴力

　　题目太长了略去不表。
　　听说这个题正解是$O(n+m\sqrt n)$的，然而我太菜只会暴力。下面来讲讲我的搞笑做法。

---

　　观察$f(a,b)$的等式，可以写成$f(a,b)=b/(a-b)f(a,b-a),b>a$，进一步的是$f(a,b)=b(a\%b)f(a,a\%b)$，若设$g=\gcd(a,b)$，则$f(a,b)=k(a,b)f(g,g)$，其中$k(a,b)$是一个和$f$值无关只与$a,b$有关的值，那么每个修改实际上也就是对$f(g,g)$的修改。
　　设$f_i=f(i,i),ans_i=\sum_{\gcd(a,b)=i}k(a,b)f_i$。考虑$f(a,b)=ab$,则$ab=k(a,b)\gcd(a,b)^2$，也即$k(a,b)=ab/\gcd(a,b)^2$。
　　因此我们有$ans_k=f_i\sum_{ik在范围内}\sum_{jk在范围内}[\gcd(i,j)==1]ij$。
　　然后。。。。。
　　如果我们可以$O(1)$获得后面的式子，那么每次修改可以$O(m)$的计算答案。。。
　　而这要求我们对$O(\min(n,m^2))$个可能的$k$计算$\sum_{i=1}^k\sum_{j=1}^k [\gcd(i,j)==1]ij$。
　　考虑递推，对于$n>1,S_n=S_{n-1}+2\sum_{i=1}^{n}[\gcd(i,n)==1]i$，设$g_n=\sum_i [\gcd(i,n)==1]i=\sum_{d\mid n}\mu(d)d(n/d+1)(n/d)/2$，这个$g_n$看起来像是积性函数和非积性函数的狄利克雷卷积，所以我们暴力一下$O(n\ln n)$地预处理，然后就可以$O(n)$算$S(n)$，然后就可以$O(m^2)$地回答$m$个询问。
　　这样的复杂度是$O(n\ln n+m^2)$，然而由于中间有大量取模/longlong乘法，所以实际效果没比暴力好多少。。。
　　我们想办法把log去掉，注意对于$n>1,g_n$实际上就是$n^2/2\sum_{d\mid n}\mu(d)/d$，这样就可以直接线性筛了。
　　复杂度就变成了$O(n+m^2)$。。。是个能A的非常优秀的暴力算法呀233333333

---

【代码】
不要问我为什么有些奇怪的无用的东西。。

#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i = a, _ = b; i <= _; i ++)
#define cerr cerr << ">>> "
const int mod = 1000000007;
const int N = 4000003;

inline int Pow(int a, int b) {
  a %= mod;
  int t = 1;
  while (b) {
    if (b & 1) t = 1ll * t * a % mod;
    a = 1ll * a * a % mod , b >>= 1;
  }
  return t;
}

inline int ars(int n) {
  return 1ll * n * (n + 1) / 2 % mod;
}

int p[N], f[N], vis[N], S[N], pr[N / 10], tot;

int e[N], mn[N], inv[N];

void init(int n) {
  inv[1] = 1;
  rep (i , 2 , n) inv[i] = mod - 1ll * (mod / i) * inv[mod % i] % mod;
  S[1] = 1;
  rep (i , 2 , n) {
    if (!vis[i])
      pr[++ tot] = i, S[i] = mod + 1 - inv[i], mn[i] = i, e[i] = i;
    for (int j = 1; j <= tot && i * pr[j] <= n; j ++) {
      int d = pr[j], t = i * d;
      vis[t] = 1;
      if (i % d == 0) {
        mn[t] = d;
        e[t] = e[i] * d; 
        S[t] = (1ll * S[t / e[t]] - 1ll * inv[d] * S[t / e[t]]) % mod;
        break;
      }
      e[t] = mn[t] = min(mn[i], d);
      S[t] = 1ll * (1 - inv[mn[t]] + mod) * S[i] % mod;
    }
  }
  int inv2 = (mod + 1) / 2;
  rep (i , 2 , n) S[i] = (1ll * S[i] * i % mod * i % mod * inv2 % mod + mod) % mod;
  S[1] = 1;
  rep (i , 2 , n) S[i] = (S[i - 1] + 2ll * i * S[i] % mod) % mod;
  rep (i , 1 , n) f[i] = 1ll * i * i % mod, p[i] = 1;
}

int x[N], g[N];
int n, m;

int gao(int a, int b) {
  if (a > b) swap(a, b);
  if (b % a == 0) return b / a % mod;
  return 1ll * b * Pow(b % a, mod - 2) % mod * gao(a, b % a) % mod;
}

int main() {
  ios::sync_with_stdio(0);
  cin >> m >> n;
  init(n);
  rep (i , 1 , m) {
    int k, a, b;
    long long y;
    cin >> a >> b >> y >> k;
    x[i]  = y % mod;
    g[i]  = __gcd(a, b);
    int tw  = gao(a, b);
    int  w  = 1ll * tw * f[g[i]] % mod;
    f[g[i]] = 1ll * x[i] * Pow(tw, mod - 2) % mod;
    int ans = 1ll * ars(k) * ars(k) % mod;
    x[i] = 1ll * x[i] * Pow(w, mod - 2) % mod;
    w = x[i];
    x[i] = 1ll * (w + mod - 1) * p[g[i]] % mod;
    p[g[i]] = 1ll * p[g[i]] * w % mod;
    rep (j , 1 , i)
      (ans += 1ll * g[j] * g[j] % mod * S[k / g[j]] % mod * x[j] % mod) %= mod;
    if (ans < 0) ans += mod;
    printf("%d\n", ans);
  }
  return 0;
}