[数论][莫比乌斯反演][杜教筛] 51Nod 1220 约数之和

$Description$

求

$$\sum_{i=1}^n\sum_{i=1}^nd(ij)$$

$Solution$

$$\begin{eqnarray}\sum_{i=1}^n\sum_{j=1}^nd(ij)&=&\sum_{i=1}^n\sum_{j=1}^n\sum_{a\mid i}\sum_{b\mid j}[(a,b)=1]{aj\over b}\\&=&\sum_{k=1}^n\mu(k)kA(\lfloor{n\over k}\rfloor)B(\lfloor{n\over k}\rfloor)\\A(n)&=&\sum_{a=1}^na\lfloor{n\over a}\rfloor\\B(n)&=&\sum_{b=1}^n\sum_{j=1}^{\lfloor{n\over b}\rfloor}j\\S(n)&=&\sum_{i=1}^n(\mu\cdot id)(i)\\&=&\sum_{i=1}^n[(\mu*1)\cdot id](i)-\sum_{i=2}^nS(\lfloor{n\over i}\rfloor)i\\&=&1-\sum_{i=2}^nS(\lfloor{n\over i}\rfloor)i\end{eqnarray}$$

$A,B$分块搞，$\mu\cdot id$杜教筛。

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int MOD = 1000000007;
const int N = 1010101;

inline char get(void) {
  static char buf[100000], *S = buf, *T = buf;
  if (S == T) {
    T = (S = buf) + fread(buf, 1, 100000, stdin);
    if (S == T) return EOF;
  }
  return *S++;
}
template<typename T>
inline void read(T &x) {
  static char c; x = 0; int sgn = 0;
  for (c = get(); c < '0' || c > '9'; c = get()) if (c == '-') sgn = 1;
  for (; c >= '0' && c <= '9'; c = get()) x = x * 10 + c - '0';
  if (sgn) x = -x;
}

int prime[N], vis[N];
int mu[N], pre[N];
int n, lim, ans, Pcnt;
map<int, int> mp;

inline int Mod(int x) {
  return (x % MOD + MOD) % MOD;
}
inline void Add(int &x, int a) {
  x += a; while (x >= MOD) x -= MOD;
}
inline int Sum1(int n) {
  return (ll)n * (n + 1) / 2 % MOD;
}
inline int Sum1(int l, int r) {
  return Mod(Sum1(r) - Sum1(l - 1));
}
inline void Pre(int n) {
  mu[1] = 1; int x;
  for (int i = 2; i <= n; i++) {
    if (!vis[i]) {
      mu[i] = -1; prime[++Pcnt] = i;
    }
    for (int j = 1; j <= Pcnt && (x = prime[j] * i) <= n; j++) {
      vis[x] = 1;
      if (i % prime[j]) mu[x] = -mu[i];
      else {
    mu[x] = 0; break;
      }
    }
  }
  for (int i = 1; i <= n; i++)
    Add(pre[i], Mod(mu[i] * i + pre[i - 1]));
}
inline int S(int n) {
  if (n <= lim) return pre[n];
  if (mp.count(n)) return mp[n];
  int res = 1, pos;
  for (int i = 2; i <= n; i = pos + 1) {
    pos = n / (n / i);
    Add(res, MOD - (ll)S(n / i) * Sum1(i, pos) % MOD);
  }
  return mp[n] = res;
}
inline int A(int n) {
  int res = 0, pos;
  for (int a = 1; a <= n; a = pos + 1) {
    pos = n / (n / a);
    Add(res, (ll)Sum1(a, pos) * (n / a) % MOD);
  }
  return res;
}
inline int B(int n) {
  int res = 0, pos;
  for(int b = 1; b <= n; b = pos + 1) {
    pos = n / (n / b);
    Add(res, (ll)Sum1(n / b) * (pos - b + 1) % MOD);
  }
  return res;
}

int main(void) {
  freopen("1.in", "r", stdin);
  freopen("1.out", "w", stdout);
  read(n); Pre(lim = pow(n, 0.666));
  int pos;
  for (int i = 1; i <= n; i = pos + 1) {
    pos = n / (n / i);
    Add(ans, (ll)Mod(S(pos) - S(i - 1)) * A(n / i) % MOD * B(n / i) % MOD);
  }
  cout << ans << endl;
  return 0;
}