This is not the end.

0%

Min_25筛

Posted on

食用方法

用于求积性函数的前$n$项的和，$n$通常在$10^{10}$左右。它要求这个函数必须是积性函数，并且对于质数，这个函数的取值可以表示成一个可以快速算前缀和并且完全积性的东西（比如低阶多项式），对于质数的k次方，这个函数大约可以在$O(1)$的时间内查询。

基本原理

以计算$[1,n]$的质数个数为例，即定义$f(x) = [x\in \text{Prime}]$，要求$\sum_{i=1}^n f(i)$。

用$p_j$表示第$j$个质数，定义$g(n,j)$表示$[2,n]$中所有的质数以及最小质因子大于$p_j$的数的函数值和（注意没有包含$1$）。我们将按照$j$从小到大的顺序依次计算$g(n,j),g(\lfloor {n\over 2 } \rfloor ,j),g(\lfloor {n\over 3} \rfloor ,j)\cdots $。这相当于是在模拟埃氏筛的过程。

最初我们有$g(n,0)=n-1$。

考虑从$g(n,j-1)$到$g(n,j)$，我们将会减去一些最小质因子为$p_j$的数的函数值。这些数必然会包含$p_j$这个因子，因此要减去的就是$f(p_j)$乘以$[2,\lfloor {n\over p_j}\rfloor ]$中最小质因子大于等于$p_j$的数的函数值的和。设$s_j$表示$\sum_{i=1}^j f(p_i)$，那么要减去的就是$f(p_j) \cdot (g(\lfloor {n\over p_j}\rfloor ,j-1) - s_{j-1})$。

考虑这一步的复杂度：首先由于我们只需要用到$\lfloor {n\over p}\rfloor $，所以我们只需要存储$\lfloor {n\over 1}\rfloor,\lfloor {n\over 2}\rfloor,\lfloor {n\over 3}\rfloor\cdots $这些点的取值，共$O(\sqrt n)$个；$j$这一维可以滚动；所以空间复杂度是$O(\sqrt n)$的。注意滚动之后，我们需要在进行计算的时候从大到小枚举$n$。时间复杂度被证明是$O({n^{3\over 4 }\over \ln n})$的。

然后我们考虑计算答案。我们设$h(n,j)$表示$[2,n]$所有最小质因子大于等于$p_j$的数（包括最小质因子大于等于$p_j$的合数以及大于等于$p_j$的质数）的$f$的取值的和，那么有：

第一个$\sum$中的限制是为了让$\lfloor{ n\over p_j^k }\rfloor >p_j$，否则$\lfloor{ n\over p_j^k }\rfloor - s_j$会出错。它等于

直接递归下去计算即可。由于几乎不会访问到相同的$(j,n)$，所以不需要记忆化（这个可以通过打表看出来）。时间复杂度据说是$O({n^{3\over 4 }\over \ln n})$？反正跑得很快就对了。

模板

loj6053 简单的函数

定义一个积性函数$f(x)$，对于质数的$k$次方满足$p^k$满足$f(p^k) =p\text{ xor }k$。求$\sum_{i=1}^n f(i) \bmod {10^9+7}$。$n\le 10^{10}$

Solution：

对于除了$2$以外的任何一个素数都是奇数，所以它们异或$1$得到自己$-1$，即对于大于$2$的素数我们有$f(p)=p-1$，套板子计算出质数的和、质数的个数就可以了。主意最后还要特殊考虑$2$的贡献。

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
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#define ll long long
using namespace std;
const int N=2e5+10;
const int mod=1e9+7;
ll n;
int lim0,lim1;
int id0[N],id1[N];
ll w[N];
int ID(ll x) { return x<=lim0?id0[x]:id1[n/x]; }
int f[N],g[N],h[N];
int sump[N];
int pri[N],num;
void getpri(int n) {
	static int s[N];
	for(int i=2;i<=n;++i) {
		if(!s[i]) pri[s[i]=++num]=i;
		for(int j=1;j<=s[i];++j)
			if(pri[j]*(ll)i<=n) s[i*pri[j]]=j;
			else break;
	}
}
int calc(ll n,int i) {
	if(pri[i]>n||n<=1) return 0;
	int p=ID(n);
	ll ans=(g[p]-h[p]-(sump[i-1]-(i-1)))%mod; if(i==1) ans=(ans+2)%mod;
	for(int j=i;j<=num&&pri[j]*(ll)pri[j]<=n;++j) {
		ll l=pri[j];
		for(int e=1;l*(ll)pri[j]<=n;++e,l*=pri[j]) {
			ans=(ans+(pri[j]^e)*(ll)calc(n/l,j+1)+(pri[j]^(e+1)))%mod;
		}
	}
	return (ans%mod+mod)%mod;
}
int main() {
	scanf("%lld",&n);
	lim0=ceil(sqrt(n));
	getpri(lim0);
	int _tot=0;
	for(ll l=1,r;l<=n;l=r+1) {
		r=n/(n/l);
		ll v=n/l;
		(v<=lim0?id0[v]:id1[n/v])=++_tot;
		w[_tot]=v;
		g[_tot]=(v%mod+2)*(ll)(v%mod-2+1)%mod*((mod+1)/2)%mod;
		h[_tot]=(v-1)%mod;
	}
	for(int i=1;i<=num;++i) {
		sump[i]=(sump[i-1]+pri[i])%mod;
		for(int j=1;j<=_tot;++j) {
			if(pri[i]*(ll)pri[i]>w[j]) break;
			int k=ID(w[j]/pri[i]);
			g[j]=(g[j]-pri[i]*(ll)(g[k]-sump[i-1])%mod)%mod;
			h[j]=(h[j]-(h[k]-(i-1)))%mod;
		}
	}
	printf("%d\n",(calc(n,1)+1)%mod);
	return 0;
}

luogu P5325

定义一个积性函数$f(x)$，对于质数的$k$次方$p^k$满足$f(p^k) = p^k(p^k -1)$，求$\sum_{i=1}^n f(i)\bmod {10^9+7}$。$n\le 10^{10}$

Solution：

对于质数这个函数满足$f(x) = x(x-1) = x^2 - x$，分成$x^2$和$x$两部分分别计算就可以了。

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
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#define ll long long
using namespace std;
const int mod=1e9+7;
int Pow(int x,int y) {
    int res=1;
    while(y) {
        if(y&1) res=res*(ll)x%mod;
        x=x*(ll)x%mod,y>>=1;
    }
    return res;
}
const int N=2e5+10;
const int inv2=(mod+1)/2,inv6=Pow(6,mod-2);
ll n;
int lim0,id0[N],id1[N],m;
ll w[N];
int g[N],h[N];
int ID(ll x) { return x<=lim0?id0[x]:id1[n/x]; }
int pri[N],num;
int sg[N],sh[N];
void getpri(int n) {
    static int s[N];
    for(int i=2;i<=n;++i) {
        if(!s[i]) pri[s[i]=++num]=i;
        for(int j=1;j<=s[i];++j)
            if(pri[j]*(ll)i<=n) s[pri[j]*i]=j;
            else break;
    }
}
inline int F(ll x) { return x%mod*(x%mod-1)%mod; }
int calc(ll n,int i) {
    if(pri[i]>n||n<=1) return 0;
    int x=ID(n),ans=((g[x]-h[x])%mod-(sg[i-1]-sh[i-1])%mod)%mod;
    for(int j=i;j<=num;++j) {
        if(pri[j]*(ll)pri[j]>n) break;
        ll l=pri[j];
        for(int e=1;l*(ll)pri[j]<=n;++e,l*=pri[j]) {
            ans=(ans+F(l)*(ll)calc(n/l,j+1)+F(l*pri[j]))%mod;
        }
    }
    return ans;
}
int main() {
    scanf("%lld",&n);
    getpri(lim0=ceil(sqrt(n)));
    for(ll l=1,r;l<=n;l=r+1) {
        ll v=n/l,z=v%mod; r=n/v;
        (v<=lim0?id0[v]:id1[n/v])=++m;
        w[m]=v;
        g[m]=(z*(z+1)%mod*(2*z+1)%mod*inv6%mod-1)%mod;
        h[m]=(z+2)*(z-2+1)%mod*inv2%mod;
    }
    for(int i=1;i<=num;++i) {
        sg[i]=(sg[i-1]+pri[i]*(ll)pri[i])%mod;
        sh[i]=(sh[i-1]+pri[i])%mod;
        for(int j=1;j<=m;++j) {
            if(pri[i]*(ll)pri[i]>w[j]) break;
            int k=ID(w[j]/pri[i]);
            g[j]=(g[j]-pri[i]*(ll)pri[i]%mod*(g[k]-sg[i-1])%mod)%mod;
            h[j]=(h[j]-pri[i]*(ll)(h[k]-sh[i-1]))%mod;
        }
    }
    printf("%d\n",((calc(n,1)+1)%mod+mod)%mod);
    return 0;
}

jzoj5594 最大真因数

求$[1,n]$的所有合数的最大真因数的和。

Solution：

一个数的最大真因数等于这个数除以它的最小质因子。

考虑直接枚举这个数的最小质因子以及这个最小质因子的幂次$p^k$，则它对答案的贡献是所有最小质因子大于$p$的数的和乘以$p^{k-1}$。这可以在min_25筛的过程中解决。

可以不用记忆化。因为我们枚举最小质因子及其幂次相当于min_25筛里面的第一层展开来。

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
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <map>
#include <vector>
#define PB push_back
#define ll long long
#define ull unsigned long long
using namespace std;
const int N=2e5+10;
int pri[N],num;
void getpri(int n) {
	static int s[N];
	for(int i=2;i<=n;++i) {
		if(!s[i]) pri[s[i]=++num]=i;
		for(int j=1;j<=s[i];++j)
			if(i*(ll)pri[j]<=n) s[i*pri[j]]=j;
			else break;
	}
}
ull g[N],sg[N],Ans;
ll w[N],n;
int m;
int id0[N],id1[N],lim;
namespace Hash {
	const int pri=1000000,N=pri+10;
	int head[N];
	vector<ll> key;
	vector<ull> val;
	vector<int> nxt;
	void init() {
		memset(head,-1,sizeof(head));
		key.clear(),val.clear(),nxt.clear();
	}
	ll id(int j,ll n) { return n*(num+1)+j; }
	ull query(ll k) {
		for(int i=head[k%pri];~i;i=nxt[i])
			if(key[i]==k) return val[i];
		return -1;
	}
	void insert(ll k,ull v) {
		key.PB(k),val.PB(v);
		nxt.PB(head[k%pri]),head[k%pri]=nxt.size()-1;
	}
	ull query(int j,ll n) { return query(id(j,n)); }
	void insert(int j,ll n,ull v) { insert(id(j,n),v); }
}		
int getid(ll x) { return x<=lim?id0[x]:id1[n/x]; }
ull calc(ll n,int j) {
	if(n<pri[j]||n<=1) return 0;
//	ull t=Hash::query(j,n); if(t!=-1) return t;
//	if(mp[j].count(n)) return mp[j][n];
	int x=getid(n); ull ans=g[x]-sg[j-1];
	for(int i=j;i<=num;++i) {
		if(pri[i]*(ll)pri[i]>n) break;
		ll l=pri[i];
		for(int e=1;l*pri[i]<=n;++e,l*=pri[i]) {
			ull tmp=calc(n/l,i+1);
			ans+=l*(tmp+pri[i]);
		//	Ans+=l/pri[i]*tmp;
		}
	}
//	Hash::insert(j,n,t);
	return ans;
	//return mp[j][n]=ans;
}	
ull solve(ll _n) {
	n=_n;
	m=0,Ans=0;
/*	for(int i=1;i<=num;++i) {
		mp[i].clear();
	}*/
	Hash::init();
	lim=ceil(sqrt(n));
	for(ll l=1,r;l<=n;l=r+1) {
		ull v=n/l; r=n/v;
		(v<=lim?id0[v]:id1[n/v])=++m;
		w[m]=v;
		if(v&1) g[m]=(v+2)*((v-1)/2);
		else g[m]=((v+2)/2)*(v-1);
	}
//	printf("%d\n",m);
//	printf("%llu\n",g[1]);
	for(int i=1;i<=num;++i) {
		sg[i]=sg[i-1]+pri[i];
		for(int j=1;j<=m;++j) {
			if(pri[i]*(ll)pri[i]>w[j]) break;
			int k=getid(w[j]/pri[i]);
			g[j]-=pri[i]*(g[k]-sg[i-1]);
		}
	}
	for(int i=1;i<=num;++i) {
		if(pri[i]*(ll)pri[i]>n) break;
		for(ll l=pri[i];l*pri[i]<=n;l*=pri[i]) {
			Ans+=l/pri[i]*calc(n/l,i+1);
			Ans+=l;
		}
	}
	return Ans;
}
int main() {
	freopen("factor.in","r",stdin);
	freopen("factor.out","w",stdout);
	getpri(100000);
	ll l,r;
	scanf("%lld%lld",&l,&r);
	printf("%llu\n",solve(r)-solve(l-1));
	return 0;
}