扯淡
min_25筛是由min_25提出的求积性函数前缀和的亚线性算法,和一个叫“扩展埃氏筛”的东西有着微妙的关系。
至于是什么关系,我也不太清楚,反正有人说很像有人说就是一个东西(雾)
这段话并不是废话
约定
为了方便后面描述,这里写一些用到的约定和符号表示,以免产生恐惧
1 1 1被开除正整数籍 也就是说“前缀和”之类的都是从 2 2 2开始,对答案所求的前缀和同样,最后手动加 1 1 1
π ( n ) \pi(n) π(n)表示 1 ∼ n 1 \sim n 1∼n中质数个数的规模(其实问题不大,后面就懂了)
p i p_i pi表示第 i i i个质数,单独的 p p p均表示质数
p r i m e prime prime表示质数集合 m i n p ( i ) minp(i) minp(i)表示 i i i的最小质因子
p c p^c pc表示一个只含一个质因子的数
流程
首先所求的函数 f ( n ) f(n) f(n)需要满足:
- 是个积性函数
- 在质数处的取值 f ( p ) f(p) f(p)是一个关于 p p p的多项式。我们把这个多项式拆成若干单项式分别计算在相加,这样就变成了" f ( p ) f(p) f(p)的值是一个关于 p p p的单项式",我们记为 f ( p ) = p k f(p)=p^k f(p)=pk
- f ( p c ) f(p^c) f(pc)的值可以快速计算
首先考虑计算这个东西
∑ i = 2 n [ i ∈ p r i m e ] i k \sum_{i=2}^n[i\in prime]i^k i=2∑n[i∈prime]ik
注意是 i k i^k ik不是 f ( i ) f(i) f(i),虽然这里还没有区别,但它们只是质数位置相等
首先既然叫扩展埃氏筛先回忆一下埃氏筛在干什么:开始写下所有正整数,然后从小到大,如果一个没被筛说明它是质数,用它筛掉后面的倍数
我们可以用这个思路计算上面的式子,先算出所有的和,再用质数筛掉所有合数项
设 g ( n , j ) g(n,j) g(n,j)表示用前 j j j个质数筛了之后剩余项的和
g ( n , j ) = ∑ i = 2 n [ i ∈ p r i m e o r m i n p ( i ) > p j ] i k g(n,j)=\sum_{i=2}^n[i \in prime \quad or\quad minp(i)>p_j]i^k g(n,j)=i=2∑n[i∈primeorminp(i)>pj]ik
和就是 g ( n , π ( n ) ) g(n,\pi(\sqrt{n})) g(n,π(n))
注意合数的最小质因子不会超过 n \sqrt n n,所以如果 p j 2 > n p_j^2>n pj2>n有 g ( n , j ) = g ( n , j − 1 ) g(n,j)=g(n,j-1) g(n,j)=g(n,j−1)
对剩下情况考虑转移,即从 g ( n , j − 1 ) g(n,j-1) g(n,j−1)去掉被 p j p_j pj删掉的项
g ( n , j − 1 ) = ∑ i = 2 n [ i ∈ p r i m e o r m i n p ( i ) ≥ p j ] g(n,j-1)=\sum_{i=2}^n[i\in prime \quad or \quad minp(i)\geq p_j] g(n,j−1)=i=2∑n[i∈primeorminp(i)≥pj]
如果脑子转不过来,可以强行算
[ i ∈ p r i m e o r m i n p ( i ) ≥ p j ] [i\in prime \quad or \quad minp(i)\geq p_j] [i∈primeorminp(i)≥pj]且不满足 [ i ∈ p r i m e o r m i n p ( i ) > p j ] [i \in prime \quad or\quad minp(i)>p_j] [i∈primeorminp(i)>pj]
即
[ i ∈ p r i m e o r m i n p ( i ) ≥ p j ] [i\in prime \quad or \quad minp(i)\geq p_j] [i∈primeorminp(i)≥pj]且 [ i ∉ p r i m e ] 且 [ m i n p ( i ) ≤ p j ] [i \notin prime] 且 [minp(i)\leq p_j] [i∈/prime]且[minp(i)≤pj]
综上
[ i ∉ p r i m e , m i n p ( i ) = p j ] [i \notin prime,minp(i)=p_j] [i∈/prime,minp(i)=pj]
所以
g ( n , j ) = g ( n , j − 1 ) − ∑ i = 2 n [ i ∉ p r i m e , m i n p ( i ) = p j ] i k g(n,j)=g(n,j-1)-\sum_{i=2}^n[i \notin prime,minp(i)=p_j]i^k g(n,j)=g(n,j−1)−i=2∑n[i∈/prime,minp(i)=pj]ik
可以提一个 p j p_j pj出来,因为不能有质数,就把 p j p_j pj去掉,刚好是 1 1 1
g ( n , j ) = g ( n , j − 1 ) − p j k ∑ i = 2 ⌊ n p j ⌋ [ m i n p ( i ) ≥ p j ] i k g(n,j)=g(n,j-1)-p_j^k\sum_{i=2}^{\lfloor\frac{n}{p_j}\rfloor}[minp(i)\geq p_j]i^k g(n,j)=g(n,j−1)−pjki=2∑⌊pjn⌋[minp(i)≥pj]ik
再次观察
g ( n , j ) = ∑ i = 2 n [ i ∈ p r i m e o r m i n p ( i ) > p j ] i k g(n,j)=\sum_{i=2}^n[i \in prime \quad or\quad minp(i)>p_j]i^k g(n,j)=i=2∑n[i∈primeorminp(i)>pj]ik
发现可以拆成小于 p j p_j pj的质数和大于等于 p j p_j pj的所有数,第二个就是刚才的式子
g ( n , j ) = g ( n , j − 1 ) − p j k ( g ( ⌊ n p j ⌋ , j − 1 ) − ∑ i = 1 j − 1 p i k ) g(n,j)=g(n,j-1)-p_j^k(g(\lfloor\frac{n}{p_j}\rfloor,j-1)-\sum_{i=1}^{j-1}p_i^k) g(n,j)=g(n,j−1)−pjk(g(⌊pjn⌋,j−1)−i=1∑j−1pik)
注意 n n n由于都是一直整除,所以只会有 O ( n ) O(\sqrt n) O(n)种取值,可以整除分块找出来强行离散化。在记录某个值 v v v所在位置的时候,如果 v > n v>\sqrt n v>n,我们另开一个数组存到 ⌊ n v ⌋ \lfloor\frac{n}{v}\rfloor ⌊vn⌋里面
然后后面只会用到最后一项,所以第二维可以滚掉
求答案时,设
S ( n , j ) = ∑ i = 2 n [ m i n p ( i ) > p j ] f ( i ) S(n,j)=\sum_{i=2}^n[minp(i)>p_j]f(i) S(n,j)=i=2∑n[minp(i)>pj]f(i)
分质数和合数分别计算
质数部分用 g g g算出来的减去前 j j j个
g ( n , π ( n ) ) − ∑ i = 1 j p i k g(n,\pi(\sqrt n))-\sum_{i=1}^jp^k_i g(n,π(n))−i=1∑jpik
合数部分枚举最小质因子和它的次数
∑ k = j + 1 p k ≤ n ∑ e = 1 p k e ≤ n f ( p k e ) ( S ( ⌊ n p k e ⌋ , k ) + [ e > 1 ] ) \sum_{k=j+1}^{p_k\leq\sqrt n}\sum_{e=1}^{p_k^e\leq n}f(p_k^e)(S(\lfloor\frac{n}{p_k^e}\rfloor,k)+[e>1]) k=j+1∑pk≤ne=1∑pke≤nf(pke)(S(⌊pken⌋,k)+[e>1])
[ e > 1 ] [e>1] [e>1]指如果指数大于 1 1 1它本身就是合数,需要统计答案
两部分相加即可
总复杂度大概 O ( n 2 3 ) O(n^{2\over3}) O(n32)常数较小,大约 2 s 2s 2s过 1 e 10 1e10 1e10, 3 s 3s 3s过 1 e 11 1e11 1e11
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cmath>
#define MAXN 200005
using namespace std;
typedef long long ll;
const int MOD=1e9+7,INV6=(MOD+1)/6;
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}
int np[MAXN],pl[MAXN],cnt;
inline void init(const int& N)
{
np[1]=1;
for (int i=2;i<=N;i++)
{
if (!np[i]) pl[++cnt]=i;
int x;
for (int j=1;(x=i*pl[j])<=N;j++)
{
np[x]=1;
if (i%pl[j]==0) break;
}
}
}
ll val[MAXN];
int tot;
int g1[MAXN],g2[MAXN],sum1[MAXN],sum2[MAXN];
int key[MAXN],yek[MAXN];
ll n;
int m;
inline int getkey(const ll& v){return v<=m? key[v]:yek[n/v];}
int S(ll n,int j)
{
if (pl[j]>=n) return 0;
int k=getkey(n);
int ans=dec(dec(g2[k],g1[k]),dec(sum2[j],sum1[j]));
for (int k=j+1;k<=cnt&&(ll)pl[k]*pl[k]<=n;k++)
for (ll e=1,pe=pl[k];pe<=n;e++,pe*=pl[k])
ans=add(ans,pe%MOD*(pe%MOD-1)%MOD*(S(n/pe,k)+(e>1))%MOD);
return ans;
}
int main()
{
scanf("%lld",&n);
m=sqrt(n);
init(m);
for (ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
val[++tot]=n/l;
int t=val[tot]%MOD;
g1[tot]=((ll)t*(t+1)/2)%MOD-1;
g2[tot]=(ll)t*(t+1)%MOD*(2*t+1)%MOD*INV6%MOD-1;
if (val[tot]<=m) key[val[tot]]=tot;
else yek[n/(n/l)]=tot;
}
for (int j=1;j<=cnt;j++)
{
for (int i=1;i<=tot&&(ll)pl[j]*pl[j]<=val[i];i++)
{
int k=getkey(val[i]/pl[j]);
g1[i]=dec(g1[i],(ll)pl[j]*dec(g1[k],sum1[j-1])%MOD);
g2[i]=dec(g2[i],(ll)pl[j]*pl[j]%MOD*dec(g2[k],sum2[j-1])%MOD);
}
sum1[j]=add(sum1[j-1],pl[j]),sum2[j]=add(sum2[j-1],(ll)pl[j]*pl[j]%MOD);
}
printf("%d\n",S(n,0)+1);
return 0;
}