Processing math: 100%

Min_25筛学习笔记

发表于 | 分类于

一直想学洲阁筛吧……（这就是我学Min_25筛的理由？）

然后今天终于能yy出来洲阁筛的复杂度证明了……（然后不想写（逃

然后发现了一种更优越的求积性函数前缀和的方法，叫做Min_25筛的……（虽然算是烂大街了？）

Min_25筛

假设现在有一个积性函数f，我们知道他在质数的幂pc上的表现f(pc)，然后要去求它的前缀和。

Min_25筛有个重要的要求就是……f(p)（p为质数）事一个低次多项式或者一个方便求前缀和的东西。下面就谈谈为什么……

我们先考虑对于所有状态⌊nx⌋，求出范围内所有质数的答案。

看起来蛮棘手的……但我们考虑模仿洲阁筛，定义一个状态g(i,x)，表示x范围内事质数或者和前i个质数都互质的数的答案之和。如果说pi达到了⌊√n⌋，那么显然g(i,x)就是x范围内所有质数的答案了（pi表示第i个质数，下同）。

那么考虑转移（接下来出现的f(pi)是指质数的答案多项式中的某一项，因为答案多项式要分开考虑）……对所有x≥p2i，使用g(i,x)=g(i−1,x)−f(pi)(g(i−1,⌊xpi⌋)−g(i−1,pi−1))转移（之所以后面又补了一个东西是因为要把多减去的含小于pi的质数的合数补回来）；至于x<p2i，我们会发现答案以后就不会变了，因为答案全部都是一堆质数的答案的和（其中不大于pi的一定不是合数，否则必定有小于等于pi的质因子；倘使大于pi的话，没有小于等于pi的质因子也只能有且仅有一个大于pi的质因子，当然就是质数了）。

这一步操作和洲阁筛几乎如出一辙，复杂度证明也完全一致（复杂度为O(n34lnn)）。

但是光考虑质数的答案没完惹……再定义状态f(i,x)表示x范围内和小于pi的质数都互质的数的答案，最终答案显然就是f(1)+f(1,n)。那么考虑一个很暴力的策略：

先把质数的答案都算进来（之前预处理了），然后考虑合数的答案。考虑枚举所有不小于pi的质数p（假设是第j个质数），如果说p2>x了那么就没法往下转移了（因为这样用上p就没法构造合数了），break出来就行了；反之则枚举p在数中所占的正指数e，对于所有pe+1≤x，对答案做f(j+1,⌊xpe⌋)f(pe)+f(pe+1)的贡献（其实就是枚举是否选够了e个p，然后考虑只用p的若干次方的情况）。

这个方法看起来极其暴力（甚至也没记忆化），但是复杂度很玄学（也很优秀）……可以参考朱老大的集训队论文。

例题：SPOJ DIVCNT3

这个题显然有f(pc)=3c+1，然后根据上面说的搞就行了……

代码：

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
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <functional>
#include <utility>
const int maxn = 400005;
using ll = long long;
int prm[maxn]; ll A[maxn], B[maxn];
ll sqr(ll x) {
  ll ret = 1;
  while((ret + 1LL) * (ret + 1LL) <= x) ret ++;
  return ret;
}

ll n, S;
int cnt;
void process() {
  S = sqr(n); cnt = 0;
#ifdef LOCAL
  printf("S : %lld\n", S);
#endif
  for(int i = 1; i <= S; i ++) {
    A[i] = i;
  }
  for(int i = 1; i <= S; i ++) {
    B[i] = n / (ll(i));
  }
  for(int i = 2; i <= S; i ++) {
    if(A[i] == A[i - 1]) continue;
    ll v = A[i - 1], lim = (ll(i)) * (ll(i));
    for(int j = 1; j <= (S / i); j ++) {
      B[j] -= B[j * i] - v;
    }
    for(int j = S / i + 1; j <= S; j ++) {
      ll th = n / (ll(j)); if(th < lim) break;
      B[j] -= A[th / (ll(i))] - v;
    }
    for(int j = S; (ll)j >= lim; j --) {
      A[j] -= A[j / (ll(i))] - v;
    }
    prm[++ cnt] = i;
  }
  prm[++ cnt] = S + 1LL;
}
inline ll query(ll x) {
  if(x <= S) {
    return A[x];
  } else {
    return B[n / x];
  }
}
ll calc(ll m, int x) {
  if(m < (ll)prm[x]) return 0;
  ll ret = 4LL * (query(m) - query(prm[x] - 1));
  for(int i = x; i <= cnt; i ++) {
    ll mul = prm[i];
    if(mul * (ll(prm[i])) > m) break;
    for(int j = 1; ; j ++) {
      if(mul * (ll(prm[i])) > m) break;
      ret += calc(m / mul, i + 1) * (ll(3 * j + 1)) + (ll(3 * j + 4));
      mul *= (ll(prm[i]));
    }
  }
  return ret;
}

int main() {
  int T; scanf("%d", &T);
  while(T --) {
    scanf("%lld", &n);
    process();
    printf("%lld\n", calc(n, 1) + 1LL);
  }
  return 0;
}