[数论][莫比乌斯反演][杜教筛] BZOJ 3512: DZY Loves Math IV

Description

求

∑i=1n∑j=1mφ(ij)

Solution

好强的题
设

S(n,m)=∑i=1mφ(ni)

当|μ(n)|=1时

φ(ni)=φ(i)∑d|(n,i)φ(nd)

这个东西为什么别人都觉得那么显然啊QAQ。。。想了一节生物课。。
令D=(i,n)，φ(i)已经把i的所有质因子的(1−1pj)全部贡献掉了，还贡献了i，后面那个东西可以考虑贡献了p∈{p是质数|(p∣n)∧(p∤i)}的(1−1pj)以及一个nD(φ∗1)(D)。。
所以就有

φ(ni)====φ(i)∑d|(n,k)φ(nd)∑i=1m∑d∣(n,k)φ(i)φ(nd)∑d∣nφ(nd)∑i=1⌊md⌋φ(di)∑d∣nφ(nd)S(d,⌊md⌋)

当μ(n)=0时令k=∏ipi那么

S(n,m)=nS(k,m)k

当n=1时就是杜教筛φ的前缀和啦。

S(1,m)=m(m+1)2−∑i=2mS(1,⌊mi⌋)

到这里就完成啦。
啊啊啊啊啊一开始我把模数和线性筛打错了QAQ

#include <bits/stdc++.h>using namespace std;typedef pair<int, int> Pairs;typedef long long ll;const int N = 1010101;const int MOD = 1000000007;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 phi[N], pre[N], mx[N], mu[N];int prime[N], vis[N];int n, m, Pcnt, lim, ans;map<Pairs, int> mp;inline void Add(int &x, int a) {  x = (x + a) % MOD;}inline void Pre(int n) {  mu[1] = 1; phi[1] = 1; mx[1] = 1; int x;  for (int i = 2; i <= n; i++) {    if (!vis[i]) {      phi[i] = i - 1; mu[i] = 1;      mx[i] = i; prime[++Pcnt] = i;    }    for (int j = 1; j <= Pcnt && (x = prime[j] * i) <= n; j++) {      vis[x] = 1;      if (i % prime[j]) {    phi[x] = phi[i] * phi[prime[j]];    mu[x] = -mu[i]; mx[x] = mx[i] * prime[j];      } else {    phi[x] = phi[i] * prime[j];    mu[x] = 0; mx[x] = mx[i]; break;      }    }  }  for (int i = 1; i <= n; i++)    pre[i] = (pre[i - 1] + phi[i]) % MOD;}inline int S(int n, int m) {  int res, pos;  if (!n || !m) return 0;  if (n == 1 && m <= lim) return pre[m];  if (mp.count(Pairs(n, m))) return mp[Pairs(n, m)];  if (n == 1) {    res = (ll)m * (m + 1) / 2 % MOD % MOD;    for (int i = 2; i <= m; i = pos + 1) {      pos = m / (m / i);      Add(res, MOD - (ll)S(1, m / i) * (pos - i + 1) % MOD);    }  } else if (mu[n] == 0){    res = (ll)n / mx[n] * S(mx[n], m) % MOD;  } else {    res = 0; int d;    for (d = 1; d * d < n; d++) {      if (n % d) continue;      Add(res, (ll)phi[n / d] * S(d, m / d) % MOD);      Add(res, (ll)phi[d] * S(n / d, m / (n / d)) % MOD);    }    if (d * d == n) Add(res, (ll)phi[d] * S(d, m / d) % MOD);  }  return mp[Pairs(n, m)] = res;}int main(void) {  freopen("1.in", "r", stdin);  freopen("1.out", "w", stdout);  read(n); read(m);  Pre(lim = max(n, (int)pow(m, 0.666)));  for (int i = 1; i <= n; i++)    Add(ans, S(i, m));  cout << ans << endl;  return 0;}