给定正整数\(n, m\)，求\(\sum_{i = 1}^n\sum_{j = 1}^m \varphi(ij)\)，膜\(10^9+7\)输出。

\(n\leq 10^5, m\leq 10^9\)。

神题……

\(n\)范围不大，似乎是在暗示我们枚举\(i\)。因此我们设状态\(S(n, m) = \sum_{i = 1}^m \varphi(in)\)。

首先有一个非常简单的结论就是：

\[\varphi(ij) = \varphi(i)\varphi(\frac{j}{\gcd(i, j)})\gcd(i, j)\]

这结论倒也显然，照着质因子搞一搞就能明白。然后下面如果反演那就凉了……

orz了标算之后发现标算用了\(\varphi\ast 1 = \mathrm{id}\)这个性质……因此有（下面假设\(n\)的所有不同质因子之积为\(w\)，\(y = \frac{n}{w}\)，至于为啥要这样你看下面柿子就懂了……）：

\[
\begin{aligned}
S(n, m)&=y\sum_{i = 1}^m \varphi(i)\varphi(\frac{w}{(i, w)})\sum_{d|(i, w)}\varphi(d)\\
&=y\sum_{i = 1}^m \varphi{i}\sum_{d|(i, w)}\varphi(\frac{w}{d})\\
&=y\sum_{d | w}\varphi(\frac{w}{d})\sum_{i = 1}^{\lfloor\frac{m}{d}\rfloor}\varphi(id)\\
&=y\sum_{d | w}\varphi(\frac{w}{d})S(d, \lfloor\frac{m}{d}\rfloor)
\end{aligned}
\]

这样一来就可以用类似于杜教筛的方法力……\(S(1, m)\)就是\(\varphi\)的前缀和，直接杜教筛，其他地方就是直接这么做就完了……粗略的估计复杂度就是\(O(m^{\frac{2}{3}} + n^{\frac{3}{4}}\sqrt{m})\)，但是注意到这个推导过程中用到了\(\sigma_0(n)\sim\sqrt{n}\)，事实上远远到不了这个级别（更何况这些\(n\)都满足每个质因子指数为\(1\)），所以跑的匪快……

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <functional>
#include <utility>
#include <vector>
#include <tr1/unordered_map>
typedef long long ll;
const ll ha = 1000000007LL;
const int N = 100000;
bool vis[N + 5]; int prm[N + 5], bs[N + 5];
ll phi[N + 5], phi_S[N + 5];
std::vector<int> dv[N + 5];
inline void process() {
  vis[1] = true; phi[1] = 1; bs[1] = 1;
  int cnt = 0;
  for(int i = 2; i <= N; i ++) {
    if(!vis[i]) {
      phi[i] = i - 1; bs[i] = i;
      prm[cnt ++] = i;
    }
    for(int j = 0; j < cnt; j ++) {
      int v = i * prm[j]; if(v > N) break;
      vis[v] = true;
      if(i % prm[j] == 0LL) {
        phi[v] = phi[i] * prm[j]; bs[v] = bs[i];
        break;
      } else {
        phi[v] = phi[i] * phi[prm[j]];
        bs[v] = bs[i] * prm[j];
      }
    }
  }
  for(int i = 1; i <= N; i ++) {
    phi_S[i] = (phi_S[i - 1] + phi[i]) % ha;
  }
  for(int i = 1; i <= N; i ++) {
    for(int j = i; j <= N; j += i) {
      dv[j].push_back(i);
    }
  }
}

const int maxn = 100005;
typedef std::pair<int, int> pii;
std::tr1::unordered_map<int, ll> d[maxn];
ll S1(ll n) {
  static const ll inv_2 = ha / 2LL + 1LL;
  ll ret = n * (n + 1LL) % ha;
  ret = ret * inv_2 % ha;
  return ret;
}
ll calc(int n, int m) {
#ifdef LOCAL
  // printf("State (%d, %d)\n", n, m);
#endif
  if(m == 0) return 0;
  if(n == 1) {
    if(m <= N) return phi_S[m];
    if(d[1].count(m)) return d[1][m];
    ll ans = S1(m);
    for(int i = 2; i <= m;) {
      int next = m / (m / i);
      ll th = (ll(next - i + 1)) * calc(1, m / i) % ha;
      ans = (ans - th + ha) % ha;
      i = next + 1;
    }
    d[1][m] = ans;
    return ans;
  }
  if(d[n].count(m)) return d[n][m];
  ll ans = 0;
  for(int i = 0; i < dv[n].size(); i ++) {
    int nn = dv[n][i];
    ll delta = phi[n / nn] * calc(nn, m / nn) % ha;
    ans = (ans + delta) % ha;
  }
  d[n][m] = ans;
  return ans;
}
ll query(int n, int m) {
  int w = bs[n]; int y = n / w;
  return (ll)y * calc(w, m) % ha;
}

int main() {
  process();
  int n, m; scanf("%d%d", &n, &m);
  ll ans = 0;
  for(int i = 1; i <= n; i ++) {
    ans = (ans + query(i, m)) % ha;
  }
  printf("%lld\n", ans);
  return 0;
}

• 相关文章
• 全局相关文章