min25 筛
0. 简介
主要解决积性函数前缀和 ( 1 0 10 ⋅ ⋅ ⋅ 1 0 11 10^{10}···10^{11} 1010⋅⋅⋅1011范围内)
积性函数定义：
{
f
(
1
)
=
1
f
(
a
b
)
=
f
(
a
)
∗
f
(
b
)
(
g
c
d
(
a
,
b
)
=
1
)
f
(
p
k
)
=
多项式
(
可快速计算
)
(
p
为素数
)
\begin{cases} f(1) = 1 \\ f(ab) = f(a) * f(b) \ \ \ (gcd(a,b) = 1) \\ f(p^k) = 多项式(可快速计算) (p为素数) \end{cases}
⎩
⎨
⎧f(1)=1f(ab)=f(a)∗f(b) (gcd(a,b)=1)f(pk)=多项式(可快速计算)(p为素数)
复杂度：时间复杂度：
O
(
n
0.75
log
n
)
O( \frac {n^{0.75}} {\log n})
O(lognn0.75) 空间复杂度：
O
(
n
)
O(\sqrt{n})
O(n)
1. 简单例子 p a r t 1 part_1 part1
为了更好的描述
m
i
n
25
min25
min25 筛的过程，我们举一个例子，以最简单的素数和为例，即：
f
(
n
)
=
∑
i
=
1
n
i
[
i
为素数
]
′
[
]
′
为逻辑符号，条件成立时为
1
，否则为
0
f(n) = \sum _{i=1} ^{n} i[i为素数] \\ '[]'为逻辑符号，条件成立时为1，否则为0
f(n)=i=1∑ni[i为素数]′[]′为逻辑符号，条件成立时为1，否则为0
m
i
n
25
min25
min25 筛借鉴了埃筛的做法，也是逐步通过素数，剔除合数。我们回顾这个过程：
- 选出一个素数
- 将这个素数的倍数全部剔除
埃筛有一个简单优化版，对于选定的素数 p r i m e prime prime，并不是从两倍开始筛除，而是从 p r i m e 2 prime^2 prime2 开始筛除，这个筛法的意义是：筛除最小质因子为 p r i m e prime prime 的合数。这个结论显然正确，同时也说明，当 p r i m e 2 > n prime^2 > n prime2>n 时，筛除过程就已经结束了。
m
i
n
25
min25
min25 筛将这个过程转化为了一个
D
P
DP
DP 过程，定义
d
p
[
n
,
j
]
dp[n,j]
dp[n,j] 表述为：
d
p
[
n
,
j
]
=
∑
i
=
1
n
i
[
i
为素数或者
i
的最小质因子大于第
j
的素数
]
dp[n,j] = \sum _{i=1} ^n i[i为素数或者i的最小质因子大于第j的素数]
dp[n,j]=i=1∑ni[i为素数或者i的最小质因子大于第j的素数]
这个
d
p
[
n
,
j
]
dp[n,j]
dp[n,j] 在埃筛中的意义为，对于一个数列，已经筛除了最小质因子小于等于
p
r
i
m
e
[
j
]
prime[j]
prime[j] 的数，即已经处理完了第
j
j
j 个质数的。在剩下的数列中只剩下素数和最小质因子大于
p
r
i
m
e
[
j
]
prime[j]
prime[j] 的合数。
例如： d p [ 25 , 2 ] dp[25,2] dp[25,2] 剩下的数列为 ( p r i m e [ 1 ] = 2 , p r i m e [ 2 ] = 3 ) (prime[1] = 2,prime[2] = 3) (prime[1]=2,prime[2]=3)
2，3，4，5，6，7，8，9，10，11，12，13，14，15，16，17，18，19，20 ，21，22，23，24，25
那么如何转移呢？显然
j
j
j 每次增加，都会筛除一部分，也就是
d
p
dp
dp 的值会减小，我们思考埃筛的转移，对于一个质数
p
r
i
m
e
prime
prime ，
1
到
n
1到n
1到n中的数
p
r
i
m
e
prime
prime 的倍数有
⌊
n
p
r
i
m
e
⌋
\lfloor \frac{n}{prime} \rfloor
⌊primen⌋ 个，换一种理解方式，可以理解为：
1
1
1 到
n
n
n 中每个的
p
r
i
m
e
prime
prime 的倍数除
p
r
i
m
e
prime
prime 后的商组成的集合，显然这个集合中，有一些数早已被小于
p
r
i
m
e
prime
prime 的素数筛除，我们还需要剔除这些数，这时我们惊奇的发现，这个剔除后的集合就是
d
p
[
⌊
n
p
r
i
m
e
⌋
,
j
−
1
]
dp[\lfloor \frac{n}{prime} \rfloor,j-1]
dp[⌊primen⌋,j−1]，但是还有一点问题，这个集合中包含了小于
p
r
i
m
e
prime
prime 的素数，这是我们不需要的，因此还需要减去
1
1
1 到
p
r
i
m
e
−
1
prime-1
prime−1 中的素数集合，即
d
p
[
p
r
i
m
e
−
1
,
j
−
1
]
dp[prime-1,j-1]
dp[prime−1,j−1] ；至此我们得到了完整的转移方程：
d
p
[
n
,
j
]
=
d
p
[
n
,
j
−
1
]
−
p
r
i
m
e
[
j
]
∗
(
d
p
[
⌊
n
p
r
i
m
e
⌋
,
j
−
1
]
−
d
p
[
p
r
i
m
e
−
1
,
j
−
1
]
)
d
p
[
n
,
0
]
=
n
∗
(
n
+
1
)
2
−
1
;
dp[n,j] = dp[n,j-1] - prime[j] * (dp[ \lfloor \frac{n}{prime} \rfloor ,j-1] - dp[prime-1,j-1]) \\ dp[n,0] = \frac {n * (n + 1)} {2} - 1;
dp[n,j]=dp[n,j−1]−prime[j]∗(dp[⌊primen⌋,j−1]−dp[prime−1,j−1])dp[n,0]=2n∗(n+1)−1;
为什么要乘以
p
r
i
m
e
[
j
]
prime[j]
prime[j] ？这是因为
(
d
p
[
⌊
n
p
r
i
m
e
⌋
,
j
−
1
]
−
d
p
[
p
r
i
m
e
−
1
,
j
−
1
]
)
(dp[ \lfloor \frac{n}{prime} \rfloor ,j-1] - dp[prime-1,j-1])
(dp[⌊primen⌋,j−1]−dp[prime−1,j−1]) 实际上是商的和，为了还原原始集合，所以乘以
p
r
i
m
e
[
j
]
prime[j]
prime[j]。
d p [ n , 0 ] dp[n,0] dp[n,0] 实际上就算 1 1 1 到 n n n 的前缀和，为什么减 1 1 1 ？这是为了在转移过程中不删除自己，也就是一倍的自己（方便了 D P DP DP 的转移）最后的结果加上 1 1 1 即可。
当然我们注意到了， n n n 的范围是 1 0 11 10^{11} 1011， j j j 的范围是 p r i m e [ j ] 2 < n prime[j]^2 < n prime[j]2<n (具体多少不知道)，显然第一维都存不下。注意到 n n n 的转移只与 ⌊ n p r i m e ⌋ \lfloor \frac{n}{prime} \rfloor ⌊primen⌋ 有关，而关于 ⌊ n x ⌋ \lfloor \frac{n}{x} \rfloor ⌊xn⌋ 的值的数量，实际上只有 n \sqrt{n} n 个，因此我们给 ⌊ n x ⌋ \lfloor \frac{n}{x} \rfloor ⌊xn⌋ 编号，第一维就压缩到了 n \sqrt{n} n 个，而第二维实际上可以推过去，类似于背包问题，不需要存。空间复杂度为 O ( n ) O(\sqrt n) O(n)。时间复杂度分析：首先第一维有一个 n \sqrt n n 的复杂度，第二维是 p r i m e j 2 < = n prime_j^2 <= n primej2<=n 的素数个数( n < 1 0 10 n < 10^{10} n<1010时不超过1w个)。时间复杂度(极限情况下1e9左右)可以承受。
其他注意事项：
我们考虑给 ⌊ n x ⌋ \lfloor \frac{n}{x} \rfloor ⌊xn⌋ 编号的时候，需要用两个数组储存编号：
void get_id() // 数论分块，求编号
{
ll i = 1,j;
while(i <= n)
{
j = n / (n / i);
w[++tot] = n / i; // w[tot] 储存的时编号为tot的值
g[tot] = w[tot] % MOD * (w[tot] % MOD + 1) / 2 % MOD; // 前缀和
g[tot]--;
if (w[tot] > Sqr) id2[n/w[tot]] = tot;
else id1[w[tot]] = tot;
/*
Sqr = sqrt(n)
两个id数组储存的值为w[tot]时的编号,以防越界，用n/w[tot]反转一下，id2储存大于sqrt的值
id1储存小于的，用于值反解编号。
*/
i = j + 1;
}
}
inline int ID(ll x) // 反解编号
{
if (x > Sqr) return id2[n/x];
else return id1[x];
}
d p dp dp的转移代码如下
void get_g()
{
for (int i = 1; (ll)prime[i] * prime[i] <= n; ++i)
{
ll temp = (ll)prime[i] * prime[i];
for (int j = 1; j <= tot && temp <= w[j]; ++j)
{
g[j] = (g[j] - prime[i] * (g[ID(w[j]/prime[i])] - g[ID(prime[i] - 1)]) + MOD) % MOD;
}
}
}
w[i]的值是从大到小的遍历的，所以最终的结果为 g[ID(n)]
这个
D
P
DP
DP 过程是非常绝妙的，但是它只是
m
i
n
25
min25
min25 筛的基本思路，这个过程主要是求的积性函数
f
(
x
)
f(x)
f(x) 在素数部分的前缀和，即：
F
(
n
)
=
∑
i
=
1
n
f
(
i
)
[
i
为素数
]
F(n) = \sum _{i=1} ^n f(i)[i为素数]
F(n)=i=1∑nf(i)[i为素数]
如果你理解了上述过程，我们将进行扩展，进行完整的
m
i
n
25
min25
min25 讲解。
2. 更进一步 p a r t 2 part_2 part2
思考一个问题，上述筛出素数的过程中定义的 d p [ n , 0 ] dp[n,0] dp[n,0] 边界，意义是什么？在素数和中定义为前缀和，实际上背后的意义是，基于对素数部分 f ( x ) f(x) f(x) 定义的另一个一个完全积性函数。可以理解为一个辅助函数 k ( x ) k(x) k(x)，这个函数是完全积性的，并且在素数部分的定义与 f ( x ) f(x) f(x) 相同。即 f ( x ) = k ( x ) , x ∈ p r i m e f(x) = k(x),x \in prime f(x)=k(x),x∈prime，我们实际上是在这个完全积性函数的基础上作埃筛，由于最后筛出的只是素数部分，并且 f ( x ) = g ( x ) f(x)=g(x) f(x)=g(x) ，所以结果是相同的。
理解了这一步，我们就知道了，这个筛法是基于一个素数部分定义相同的完全积性函数，我们便引出了辅助函数的概念，当然这个对于后面的步骤没有太大意义，只是把原理弄清楚，方便理解。
来看一个几个变种，加深理解：
-
素数平方和，即 F ( n ) = ∑ i = 1 n i 2 [ i ∈ p r i m e ] F(n) = \sum _{i=1} ^n i^2[i\in prime] F(n)=∑i=1ni2[i∈prime] ，我们可以定义完全积性函数 f ( x ) = x 2 f(x) = x^2 f(x)=x2，按照上述 d p dp dp 转移即可，定义 d p [ n , 0 ] = n ( n + 1 ) ( 2 n + 1 ) 6 − 1 dp[n,0] = \frac {n (n+1) (2n+1)} {6} - 1 dp[n,0]=6n(n+1)(2n+1)−1。
-
素数个数，即 F ( n ) = ∑ i = 1 n 1 [ i ∈ p r i m e ] F(n) = \sum _{i=1} ^n 1 [i \in prime] F(n)=∑i=1n1[i∈prime]，我们定义完全积性函数 f ( x ) = 1 f(x) = 1 f(x)=1 ，定义 d p [ n , 0 ] = n − 1 dp[n,0] = n-1 dp[n,0]=n−1。
这个过程理解了之后，我们来考虑最初的问题，求完整的积性前缀和 (不只有素数部分)，即：
F
(
n
)
=
∑
i
=
1
n
f
(
i
)
f
(
x
)
为积性函数
F(n) = \sum _{i=1} ^n f(i) \\ f(x) 为积性函数
F(n)=i=1∑nf(i)f(x)为积性函数
积性函数具有的性质 ： $ f(ab) = f(a) * f(b),gcd(a,b) = 1成立时$ 。那么我们可以利用积性函数的性质，用素数部分前缀和还原到整个区域。
我们定义 S ( n , j ) = ∑ i = 1 n f ( i ) [ i 的最小质因子大于 p i r m e j ] S(n,j) = \sum _{i=1} ^n f(i) [i的最小质因子大于pirme_j] S(n,j)=∑i=1nf(i)[i的最小质因子大于pirmej] ，文字表述为： 当 i i i 的最小质因子大于第 j j j 个质数，且 $i <= n $ 时，计算 f ( i ) f(i) f(i) 的贡献到 S ( n , j ) S(n,j) S(n,j) 中，注意，这里没有额外包括素数，也就是说第 1 1 1 到 j j j 个质数不计入贡献。
当然我们不能忘了积性函数的另一个条件： f ( p k ) = 关于 p 的多项式 ( 可快速计算 ) ( p 为素数 ) f(p^k) = 关于p的多项式(可快速计算) (p为素数) f(pk)=关于p的多项式(可快速计算)(p为素数)， f ( p k ) f(p^k) f(pk) 为一个多项式，需要单独处理。
那么我们该如何转移呢？
我们可以枚举
j
j
j 之后的质数，并枚举质数的系数，
S
(
n
,
j
)
S(n,j)
S(n,j) 的定义为我们提供了互质的条件，因此有如下转移：
S
(
n
,
j
)
=
∑
i
=
j
+
1
∞
∑
e
=
1
∞
S
(
⌊
n
p
r
i
m
e
i
e
⌋
,
i
)
∗
f
(
p
r
i
m
e
i
e
)
S(n,j) = \sum _{i=j+1} ^{\infty} \sum _{e = 1} ^{\infty} S( \lfloor \frac{n}{prime_i^e} \rfloor,i) * f(prime_i^e)
S(n,j)=i=j+1∑∞e=1∑∞S(⌊primeien⌋,i)∗f(primeie)
当然这个转移方程并不是真正的无穷大，显然
p
r
i
m
e
i
2
>
n
prime_i^2 > n
primei2>n 时可以停止，同时
n
<
p
r
i
m
e
i
e
n < prime_i^e
n<primeie 时也可以停止。为什么这个是正确的。我们分析一下：
- S ( ⌊ n p r i m e i e ⌋ , i ) S( \lfloor \frac{n}{prime_i^e} \rfloor,i) S(⌊primeien⌋,i) 保证在 $1 到 \lfloor \frac{n}{prime_i^e} \rfloor $ 范围内所有数的最小质因子大于 p r i m e i prime_i primei ，也是变相说明了 S ( ⌊ n p r i m e i e ⌋ , i ) S( \lfloor \frac{n}{prime_i^e} \rfloor,i) S(⌊primeien⌋,i) 包含的数与 p r i m e i e prime_i^e primeie 互质，将 S 乘以 f ( p r i m e i e ) S乘以f(prime_i^e) S乘以f(primeie) ，说明在积性函数的定义下，还原了以 p r i m e i prime_i primei 为最小质因子的所有数
- 在 S S S 的定义中，并不包含 1 1 1 ，也就说，这个转移式并没有包含 f ( p r i m e i e ) f(prime_i^e) f(primeie) 这单独一项，所以我们需要加上。
具体代码实现， S ( n , j ) S(n,j) S(n,j) 我们采用递归实现：
我们以积性函数 f ( p k ) = p k ( p k − 1 ) f(p^k) = p^k(p^k-1) f(pk)=pk(pk−1) 为例子，我们可以展开为函数 f ( x ) = x 2 − x f(x) = x^2 - x f(x)=x2−x。分开处理 g 2 ( n ) = ∑ i = 1 n i 2 g2(n) = \sum _{i=1} ^n i^2 g2(n)=∑i=1ni2， g 1 ( n ) = ∑ i = 1 n i g1(n) = \sum _{i=1} ^n i g1(n)=∑i=1ni， f ( x ) = g 2 ( x ) − g 1 ( x ) f(x) = g2(x) - g1(x) f(x)=g2(x)−g1(x)，先处理出 g 1 ( x ) , g 2 ( x ) g1(x),g2(x) g1(x),g2(x) 在素数部分的和，再合并求完整积性函数。
ll S(ll x,int i) // 由素数还原整个范围
{
if (prime[i] >= x) return 0;
int id = ID(x);
ll ans = (g2[id] - g1[id] - (g2[ID(prime[i])] - g1[ID(prime[i])])) % MOD + MOD;
ans %= MOD; // 先得出大于prime[i]的素数函数和
for (int k = i+1; k <= op && (ll)prime[k] * prime[k] <= x; ++k)
// prime[k]^2>x时，S()显然等于0
{
ll pe = prime[k];
for (int e = 1; pe <= x; e++,pe *= (ll)prime[k]) // 枚举系数
{
ll t = pe % MOD;
ans += t * (t - 1) % MOD * (S(x/pe,k) + (e != 1)) % MOD;
/*
e==1时为素数，一开始就已经算入贡献了，
e > 1时，因为S() 中没有1，所以需要单独算。
注：f(t) = t * (t - 1)
*/
ans %= MOD;
}
}
return ans;
}
3. 实例讲解 p a r t 3 part_3 part3
为了更好的理解，我们给出一个实例，luogu5325模板，积性函数定义如下：
{
f
(
1
)
=
1
f
(
p
k
)
=
p
k
(
p
k
−
1
)
\begin{cases} f(1) = 1 \\ f(p^k) = p^k(p^k-1) \end{cases}
{f(1)=1f(pk)=pk(pk−1)
求：
F
(
n
)
=
∑
i
=
1
n
f
(
x
)
F(n) = \sum _{i=1} ^n f(x)
F(n)=∑i=1nf(x)
处理问题为一下三步：
- 构造素数部分完全相同的完全积性函数
- d p dp dp 转移筛选出素数部分和
- 由互质关系还原到整个定义域
我们考虑第一步，构建完全积性函数，我们可以展开为函数 f ( x ) = x 2 − x f(x) = x^2 - x f(x)=x2−x。分开处理 g 2 ( n ) = ∑ i = 1 n i 2 g2(n) = \sum _{i=1} ^n i^2 g2(n)=∑i=1ni2， g 1 ( n ) = ∑ i = 1 n i g1(n) = \sum _{i=1} ^n i g1(n)=∑i=1ni，这样我们就构建出了两个完全积性函数。
然后是第二步，我们完全可以先处理出 g 1 ( x ) , g 2 ( x ) g1(x),g2(x) g1(x),g2(x) 在素数部分的和，再合并 f ( x ) = g 2 ( x ) − g 1 ( x ) f(x) = g2(x) - g1(x) f(x)=g2(x)−g1(x)。我们就求出了 f ( x ) f(x) f(x) 在素数部分的和。
第三步，按照 f ( x ) f(x) f(x) 的定义转移即可。
贴出代码：
#include <cstdio>
#include <cmath>
using namespace std;
typedef long long ll;
const ll lim = 2e5+3;
const ll MOD = 1e9+7;
const ll inv3 = 333333336;
ll g1[lim],g2[lim],w[lim],n,Sqr;
int id1[lim],id2[lim],tot;
int prime[lim],op;
bool np[lim];
void get_prime() // 打素数表，处理素数部分函数的和
{
Sqr = sqrt(n);
for (ll i = 2; i <= Sqr; ++i)
{
if (!np[i]) prime[++op] = i;
for (ll j = 1; j <= op && prime[j] * i <= Sqr; ++j)
{
np[i * prime[j]] = 1;
if (i % prime[j] == 0)
break;
}
}
// printf("%d\n",op);
}
void get_id() // 数论分块，完全积性函数的情况下求全部的和
{
ll i = 1,j;
while(i <= n)
{
j = n / (n / i);
w[++tot] = n / i;
g1[tot] = w[tot] % MOD * (w[tot] % MOD + 1) / 2 % MOD;
g2[tot] = g1[tot] * (w[tot] * 2 % MOD + 1) % MOD * inv3 % MOD;
g1[tot]--;
g2[tot]--;
if (w[tot] > Sqr) id2[j] = tot;
else id1[w[tot]] = tot;
i = j + 1;
}
// printf("%d\n",tot);
}
inline ll ID(ll x)
{
if (x > Sqr) return id2[n/x];
return id1[x];
}
void get_g() // 由完全积性函数推出整个素数范围
{
for (ll i = 1; i <= op; ++i)
{
for (ll j = 1; j <= tot && (ll)prime[i] * prime[i] <= w[j]; ++j)
{
ll k = ID(w[j] / prime[i]);
g1[j] = (g1[j] - (ll)prime[i] * (g1[k] - g1[ID(prime[i]-1)] + MOD) % MOD + MOD) % MOD;
g2[j] = (g2[j] - (ll)prime[i] * prime[i] % MOD * (g2[k] - g2[ID(prime[i]-1)] + MOD) % MOD + MOD) % MOD;
}
}
}
ll S(ll x,int i) // 由素数还原整个范围
{
if (prime[i] >= x) return 0;
int id = ID(x);
ll ans = (g2[id] - g1[id] - (g2[ID(prime[i])] - g1[ID(prime[i])])) % MOD + MOD;
ans %= MOD;
for (int k = i+1; k <= op && (ll)prime[k] * prime[k] <= x; ++k)
{
ll pe = prime[k];
for (int e = 1; pe <= x; e++,pe *= (ll)prime[k])
{
ll t = pe % MOD;
ans += t * (t - 1) % MOD * (S(x/pe,k) + (e != 1)) % MOD;
ans %= MOD;
}
}
return ans;
}
int main()
{
scanf("%lld",&n);
get_prime();
get_id();
get_g();
printf("%lld\n",(S(n,0) + 1) % MOD);
return 0;
}
4. 更难的应用 p a r t 4 part_4 part4
定义 H a s s e Hasse Hasse 图，给定一个标准值 n n n ，定义为一个点权集合 { x ∣ x ∈ n 的因子 } \{x|x \in n的因子\} {x∣x∈n的因子}，当两个点存在连边，当且仅当两个点之间存在整除关系，并且为素数倍。 f ( n ) f(n) f(n) 定义为标准值为 n n n 的 H a s s e Hasse Hasse 的边数。要求 F ( n ) = ∑ i = 1 n f ( i ) F(n) = \sum _{i=1} ^n f(i) F(n)=∑i=1nf(i)， n < = 1 0 10 n <= 10^{10} n<=1010。
这个问题比较难了，来源于 “2021牛客多校6-G题” ，赛时共37队过题。
这个问题，在形式上比较迷惑，唯一能跟
m
i
n
25
min25
min25 靠边的只有前缀和这个条件。为了更好的分析问题，我们先分析特殊情况，比较显然的结论：
f
(
n
)
=
{
0
,
n
=
1
e
,
n
=
p
e
(
p
∈
p
r
i
m
e
)
f(n) = \begin{cases} 0 , \quad n = 1 \\ e , \quad n = p^e(p \in prime) \end{cases}
f(n)={0,n=1e,n=pe(p∈prime)
我们便定义出了在 p e p^e pe 和 1 1 1 处的函数值， m i n 25 min25 min25 显然可以处理出这个部分，那么其他部分呢？假设我们已经得到了一个任意数的值 f ( a ) f(a) f(a) ，回顾 m i n 25 min25 min25 ，由素数还原到整个域的过程，我们是枚举一个质数 p e p^e pe ,有 g c d ( p e , a ) = 1 gcd(p^e,a) = 1 gcd(pe,a)=1 ，我们便可以通过积性函数得到 f ( a ∗ p e ) f(a*p^e) f(a∗pe) ，如果我们可以得到这个转移式，这个问题即解决了。
为了得到转移式子，我们来图形化分析这个问题，我们把已经得到的 f ( a ) f(a) f(a) 想象为平面上的一幅图，我们把这个图上每个节点都乘以一个 p i , i ∈ [ 0 , e ] p^i,i \in [0,e] pi,i∈[0,e] ，显然可以得到 e + 1 e+1 e+1 幅图，每幅图都是合法的，那么至少有 ( e + 1 ) ∗ f ( a ) (e+1)*f(a) (e+1)∗f(a) 是应该被算入贡献，
第二部分，我们定义 d ( a ) d(a) d(a) 为 a a a 的因子个数，即图的节点个数，显然每幅图之间还存在节点对应连边，举个例子： x x x为 a a a的一个因子， x ∗ p i − 1 x*p^{i-1} x∗pi−1 与 x ∗ p i x*p^i x∗pi 显然存在一条连边，且存在于两幅图中。（可以理解为分层图。相邻层之间又有连边。）
这样我们就得到了一个转移式：
f
(
a
∗
p
e
)
=
(
e
+
1
)
∗
f
(
a
)
+
e
∗
d
(
a
)
f(a*p^e) = (e+1)* f(a) + e * d(a)
f(a∗pe)=(e+1)∗f(a)+e∗d(a)
这个形式虽然不是是积性函数，但是却可以用
m
i
n
25
min25
min25 解决，因为在互质的数中存在一个函数转移，并且可以得证：
s
=
a
b
,
g
c
d
(
a
,
b
)
=
1
s=ab,gcd(a,b)=1
s=ab,gcd(a,b)=1，
a
,
b
a,b
a,b不同但是
s
s
s 相同时，转移得到的
f
(
s
)
f(s)
f(s) 是相同的（从图的角度考虑证明），那么就可以用
m
i
n
25
min25
min25 解决，筛选出素数部分（显然前缀和就是素数个数），按照定义还原即可。且
d
(
n
)
d(n)
d(n) 函数（因子函数）又是一个积性函数，两者可同步计算。
贴上代码（付费题目，可以考虑与下列代码对拍验证结果的正确性）：
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
typedef long long ll;
const int lim = 2e5+3;
const int MOD = 1145140019;
ll w[lim],g[lim],Sqr,n;
int tot,op,prime[lim],id1[lim],id2[lim];
bool np[lim];
#define x first
#define y second
typedef pair<ll,ll> PR;
void get_prime()
{
for (int i = 2; i < lim; ++i)
{
if (!np[i]) prime[++op] = i;
for (int j = 1; j <= op && i * prime[j] < lim; ++j)
{
np[i * prime[j]] = 1;
if (i % prime[j] == 0)
break;
}
}
}
inline int ID(ll x)
{
if (x > Sqr) return id2[n/x];
else return id1[x];
}
void get_id()
{
Sqr = sqrt(n);
tot = 0;
ll i = 1, j;
while(i <= n)
{
j = n / (n / i);
w[++tot] = n/i;
g[tot] = (w[tot] - 1) % MOD;
if (w[tot] > Sqr) id2[j] = tot;
else id1[w[tot]] = tot;
i = j + 1;
}
}
void get_g()
{
for (int i = 1; (ll)prime[i] * prime[i] <= n; ++i)
{
ll temp = (ll)prime[i] * prime[i];
for (int j = 1; j <= tot && temp <= w[j]; ++j)
{
g[j] = (g[j] - g[ID(w[j]/prime[i])] + g[ID(prime[i] - 1)] + MOD) % MOD;
}
}
}
PR solve(ll now,int i) // 同时计算两个值，可加速，避免超时
{
PR res = {0ll,0ll},t;
if (prime[i] >= now) return res;
int k = ID(now);
/* 两种写法1
if ((ll)prime[i+1] * prime[i+1] > now)
{
res.x = g[k] - g[ID(prime[i])];
res.y = res.x * 2 % MOD;
}
else
{
res = solve(now,i+1);
ll p = prime[i+1];
for (int e = 1; p <= now; p*=prime[i+1],e++)
{
t = solve(now/p,i+1);
res.x = res.x + e + (e + 1) * t.x % MOD + e * t.y % MOD;
res.y = e + 1 + res.y + (e + 1) * t.y % MOD;
while(res.x >= MOD) res.x -= MOD;
while(res.y >= MOD) res.y -= MOD;
}
}
*/
res.x = g[k] - g[ID(prime[i])];
res.y = res.x * 2 % MOD;
for (int j = i + 1; (ll)prime[j] * prime[j] <= now; ++j)
{
ll p = prime[j];
for (int e = 1; p <= now; p*=prime[j],e++)
{
t = solve(now/p,j);
res.x = res.x + (e!=1?e:0) + (e + 1) * t.x % MOD + e * t.y % MOD;
res.y = (e!=1?e+1:0) + res.y + (e + 1) * t.y % MOD;
while(res.x >= MOD) res.x -= MOD;
while(res.y >= MOD) res.y -= MOD;
}
}
return res;
}
int main()
{
int T;
scanf("%d",&T);
get_prime();
while(T--)
{
scanf("%lld",&n);
get_id();
get_g();
printf("%lld\n",solve(n,0).x);
}
return 0;
}
5. 总结
m i n 25 min25 min25 是用于解决积性函数求前缀和问题，分三步：
- 观察积性函数形式（可能需要变换）
- 构建完全积性函数，求出素数部分前缀和
- 利用积性由素数部分还原到整个域
对于类积性函数也可求解：
f
(
n
)
=
{
0
,
n
=
1
多项式
,
n
=
p
e
(
p
∈
p
r
i
m
e
)
关于
a
,
f
(
b
)
的多项式
,
n
=
a
b
,
g
c
d
(
a
,
b
)
=
1
f(n) = \begin{cases} 0 , \quad n = 1 \\ 多项式 , \quad n = p^e(p \in prime) \\ 关于a,f(b)的多项式, \quad n = ab,gcd(a,b)=1 \end{cases}
f(n)=⎩
⎨
⎧0,n=1多项式,n=pe(p∈prime)关于a,f(b)的多项式,n=ab,gcd(a,b)=1
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
