Min-25筛学习笔记

置顶 | Posted on | Edited on | In | Views:

前言

在yx大佬不厌其烦的帮助下终于还是学会了Min-25筛。
其实不算特别难理解，主要是原先没人给讲解自己瞎走了许多弯路，代码实现也相对容易，听说可以代替杜教筛来求积性函数前缀和，名字好奇怪啊，不知道怎么来的
复杂度
求单点时略优于杜教筛。
update: Min-25是一位日本选手，精通数论，目前SPOJ rank10 orz，感谢daniel14311531

Min-25筛

第一部分

首先，如果要用Min-25筛有前置条件，即积性函数$f(i)$在i为质数时是个低次多项式，并且在i是质数的幂次时可以快速求得。
既然$f(i)$在i是质数时是个低次多项式，我们不妨把这个多项式拆开，形如
显然每个幂次不会互相干扰，那么我们就分开处理。
假设现在$f(i)=i^k(i\in prime)$
并且$f$为完全积性函数，因为合数会被筛掉，所以没有影响。
定义$P_i$表示第$i$小的质数，$sum_i$表示$\sum\limits_{j=1}^{i}f(P_j)$，$Minn_i$表示$i$的最小质因数
$g(n,i)表示\sum\limits_{j=2}^{n}[Minn_j>P_i或j\in prime]f(j)$
可以理解为用前$i$个质数做埃氏筛后剩下的数的$f$的和。
考虑如何计算$g(n,i)$。
由于是个埃氏筛的过程，我们考虑第i个质数会筛掉哪些合数。
显然会筛掉最小的合数为$P_i^2$
如果$P_i^2>n$，那么啥都筛不掉，$g(n,i)=g(n,i-1)$
如果$P_i^2\le n$，那么$P_i的[P_i,\frac{n}{P_i}]$倍都会被筛掉。
$[P_i,\frac{n}{P_i}]$这些答案存在$g(\frac{n}{P_i},i-1)$中，但是这里面还包括了$P_1,P_2…P_{i-1}$这些质数。
所以$g(n,i)=g(n,i-1)-f(P_i)*(g(\frac{n}{P_i},i-1)-sum_{i-1})$
发现n的取值都形如$\lfloor \frac{a}{b}\rfloor$，通过整数分块，预处理这$2*\sqrt{n}$个值的g即可。
利用滚动数组完成递推

第二部分

我们令$S(n,m)$表示$\sum\limits_{i=2}^{n}[Minn_i\ge P_m]f(j)$
我们要求的$S(n,1)$
考虑质数的贡献，显然是$g(n,\infty)-sum_{m-1}$
考虑合数的贡献，不能从$g$中直接得到，因为$g$中对于合数我们假设他为完全积性函数，而事实上并不是，这样做只是方便计算质数。
我们枚举每个合数的最小质因数，已及最小质因数的次数(否则无法积性函数)，那么贡献为

结合定义应该可以理解，与求g的过程相似。
S递归解决即可。

例题

Luogu 4213 Sum
让你求和
我们以求$\phi$为例，$f(i)=i-1(i\in prime)$，我们把按照上面的说法分成$f(i)=i和f(i)=-1$分别把g求出来，合并后求S即可。
求$\mu$没有本质的区别。

代码

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
// luogu-judger-enable-o2
#include<cstdio>
#include<algorithm>
#include<cctype>
#include<cstring>
#include<iostream>
#include<cmath>
#define LL long long
#define N (1000005)
using namespace std;
int T,n,pri,t,unit;
int prime[N],q[N],id1[N],id2[N];
LL g[N],h[N],sum[N];
bool nprime[N]; 
template <typename T> void read(T&t) {
    t=0;
    bool fl=true;
    char p=getchar();
    while (!isdigit(p)) {
        if (p=='-') fl=false;
        p=getchar();
    }
    do {
        (t*=10)+=p-48;p=getchar();
    }while (isdigit(p));
    if (!fl) t=-t;
}
inline void init(){
    pri=t=0;
    unit=sqrt(n);
    for (int i=2;i<=unit;i++){
        if (!nprime[i]) prime[++pri]=i,sum[pri]=sum[pri-1]+i;
        for (int j=1;j<=pri&&i*prime[j]<=unit;j++){
            nprime[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    }
    for (int i=1;i<=n;){
        int v=n/i,R=n/v;
        q[++t]=v;
        if (v<=unit) id1[v]=t;
        else id2[R]=t;
        g[t]=1ll*v*(v+1)/2-1;
        h[t]=v-1;
        i=R+1;
    }
}
inline int id(int x){
    return (x<=unit)?id1[x]:id2[n/x];
} 
LL S(int n,int m){
    if (n<=1||prime[m]>n) return 0;
    LL ret=g[id(n)]-sum[m-1]+m-1;
    for (int k=m;prime[k]*prime[k]<=n&&k<=pri;k++){
        for (int v=prime[k],p1=prime[k]-1;1ll*v*prime[k]<=n;v=v*prime[k],p1=p1*prime[k])
        ret+=(S(n/v,k+1)+prime[k])*p1; 
    }
    return ret;
}
LL D(int n,int m){
    if (n<=1||prime[m]>n) return 0;
    LL ret=h[id(n)]+(m-1);
    for (int k=m;prime[k]*prime[k]<=n&&k<=pri;k++){
        ret-=D(n/prime[k],k+1);
    } 
    return ret;
}
signed main(){
    read(T);
    while (T--){
        read(n);
        init();
        for (int i=1;i<=pri;i++){
            int v=prime[i]*prime[i];
            for (int j=1;j<=t&&v<=q[j];j++){
                int k=id(q[j]/prime[i]);
                g[j]-=1ll*prime[i]*(g[k]-sum[i-1]);
                h[j]-=h[k]-i+1;
            }
        }
        for (int i=1;i<=t;i++) g[i]-=h[i],h[i]=-h[i];
        printf("%lld %lld\n",S(n,1)+1,D(n,1)+1);
    }
    return 0;
}