Daniel_yuan's Blog

芙卡米天下第一可爱 n(*≧▽≦*)n

0%

「LibreOJ Round #11」Misaka Network 与求和 题解

发表于
本文字数： 3.9k 阅读时长 ≈ 4 分钟

求 \(\sum_{i=1}^N\sum_{j=1}^Nf^k(\gcd(i,j))\)，其中 \(f(x)\) 表示 \(x\) 的第二大质因子，每个质因子算多次，即 \(f(4)=f(6)=2\)，定义 \(f(1)=0,f(Prime)=1\)。

众所周知 Min_25 筛的本质是用 DP 求出所有 \(S(\lfloor\frac{n}{1}\rfloor),S(\lfloor\frac{n}{2}\rfloor)...\) 这 \(2\sqrt{n}\) 个值。如果可以较好的处理每个数质因子之间的关系的话，是不要求所筛函数是积性函数的，这个题就是一个典型例子。

直接设 \(f^k(n)=\sum_{d|n}g(d)\)，莫比乌斯反演得到 \(g(n)=\sum_{d|n}f^k(d)\mu(\frac{n}{d})\)。原式变成 \(\sum_{d=1}^ng(d)\lfloor\frac{n}{d}\rfloor^2\)。

后面直接整除分块，现在就只需要求 \(g\) 的前缀和。

设 \(S(n)=\sum_{i=1}^ng(i)\)，那么 \(S(n)=\sum_{i=1}^n\sum_{d|n}\mu(d)f^k(\frac{n}{d})=\sum_{d=1}^n\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f^k(i)\)。后面的求和式又可以整除分块。而前面求 \(\mu\) 的前缀和可以用杜教筛（应该也只能用杜教筛）。

考虑怎么求 \(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f^k(i)\)。可以发现 Min_25 筛求积性函数和的时候，每次是剥离了一个数的最小质因子。在这个题中，如果这个最小质因子恰好是次大质因子，那么就可以直接计算贡献，否则贡献就已经在之前算过了。所以我们设 \(F(n,x)\) 表示 \(n\) 以内所有质因子大于等于 \(Prime_x\) 的合数的 \(f^k(x)\) 的和，\(G(n)\) 表示 \(n\) 以内的质数个数。类似于 Min_25 筛的，\(F,G\) 只有 \(2\sqrt{n}\) 个值。

\(G\) 是一个经典的 DP，在此就不再赘述。

考虑求 \(F\)，有 \(F(n,x)=F(n,x+1)+\sum_{e=1}\left(F(\frac{n}{p_x^e},x+1)+p_x^k\left(G(\frac{n}{p_x^e})-(x-1)\right)\right)\)，即分最小质因子大于 \(x\) 和最小质因子恰好为 \(x\) 分别处理。

最后的真实的 \(\sum_{i=1}^nf^k(i)\) 就是 \(F(n)+G(n)\)。

总复杂度暂时还不会分析，但是 Min_25 的预处理是约 \(O(n^\frac{3}{4})\)，真实求值的时候要用两个整除分块加杜教筛，复杂度可能会达到 \(O(n^{\frac{7}{8}})\)（当然这是瞎扯）。总之能过。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
#include <bits/stdc++.h>
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define RI register int
typedef long long LL;

#define FILEIO(name) freopen(name".in", "r", stdin), freopen(name".out", "w", stdout);

using namespace std;

int const MAXN = 1e6 + 5;
int n, k, sq;
int Prime[MAXN], tot, miu[MAXN], Smiu[MAXN];
char Notprime[MAXN];

void Euler(int Max) {
  miu[1] = Smiu[1] = 1;
  for (RI i = 2; i <= Max; ++i) {
    if (!Notprime[i])
      Prime[++tot] = i, miu[i] = -1;
    for (RI j = 1; j <= tot && i * Prime[j] <= Max; ++j) {
      Notprime[i * Prime[j]] = 1;
      if (i % Prime[j] == 0) break;
      miu[i * Prime[j]] = miu[i] * -1;
    }
    Smiu[i] = Smiu[i - 1] + miu[i];
  }
}
namespace DuDuDu {
  map <int, unsigned int> mp;
  unsigned int S(int n) {
    if (n <= MAXN - 5) return Smiu[n];
    if (mp.count(n)) return mp[n];
    unsigned ans = 1;
    for (RI l = 2, r; l <= n; l = r + 1) {
      r = n / (n / l);
      ans -= (r - l + 1) * S(n / l);
    }
    return mp[n] = ans;
  }
}
namespace InitF {
// g(x) = 1,            G(n, x) = \sum g(i)  i is prime || minprime[i] > x
// f(x) = second prime, F(n, x) = \sum f(i)  i is not prime && minprime[i] >= x

  int const MX = 1e5 + 5;
  int id1[MX], id2[MX], val[MX];
  unsigned int G[MX], F[MX], Ans[MX];
  inline int Getid(int x) { return x <= sq ? id1[x] : id2[n / x]; }
  unsigned int qpow(unsigned int a, int k) {
    unsigned int re = 1;
    for (; k; k >>= 1, a *= a)
      if (k & 1) re *= a;
    return re;
  }
  void Work() {
    int cnt = 0;
    for (RI l = 1, r; l <= n; l = r + 1) {
      r = n / (n / l);
      int x = (n / l);
      if (x <= sq) id1[x] = ++cnt, val[cnt] = x;
      else     id2[n / x] = ++cnt, val[cnt] = x;
    }
    for (RI i = 1; i <= cnt; ++i)
      G[i] = val[i] - 1;
    int x;
    for (x = 1; Prime[x] * Prime[x] <= val[1]; ++x)
      for (RI i = 1; i <= cnt; ++i) {
        if (Prime[x] * Prime[x] > val[i]) break;
        G[i] = G[i] - (G[Getid(val[i] / Prime[x])] - (x - 1)); // (x - 1) = SumG[x - 1]
      }
    for (; x; --x) {
      for (RI i = 1; i <= cnt; ++i) {
        LL cur = Prime[x];
        while (cur * Prime[x] <= val[i]) {
          F[i] += F[Getid(val[i] / cur)];
          F[i] += (G[Getid(val[i] / cur)] - (x - 1)) * qpow(Prime[x], k); // (x - 1) = SumG[x - 1]
          cur *= Prime[x];
        }
      }
    }
    for (RI i = 1; i <= cnt; ++i)
      Ans[i] = F[i] + G[i];
  }
}
unsigned int SG(int n) {
  unsigned int ans = 0;
  for (RI l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    unsigned int L = DuDuDu :: S(l - 1), R = DuDuDu :: S(r);
    ans += (R - L) * (InitF :: Ans[InitF :: Getid(n / l)]);
  }
  return ans;
}

int main() {
  
#ifdef LOCAL
  FILEIO("a");
#endif

  Euler(1e6);
  cin >> n >> k;
  sq = sqrt(n);
  InitF :: Work();
  unsigned int ans = 0;
  for (RI l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    unsigned int L = SG(l - 1), R = SG(r);
    ans += (R - L) * (n / l) * (n / l);
  }
  cout << ans << endl;

  return 0;
}

// created by Daniel yuan
/*
     ________
    /        \
   / /      \ \
  / /        \ \
  \            /
   \  ______  /
    \________/
*/