My Blogs
杜教筛是一种能在 \(\mathcal O(n^{\frac 2 3})\) 的时间复杂度内求积性函数前缀和的筛法。虽然复杂度比较优秀，但是被筛的积性函数需要满足特殊性质。
Min_25 筛由 Min_25 发明，相对更通用，其时间复杂度为 \(\mathcal O(\frac{n^{\frac 3 4}}{\log n})\)。
首先构造一个完全积性函数 \(f'(x)\) 满足其在质数处和 \(f(x)\) 取值相同，其余位置任意。下设 \(\mathbb{P}\) 为质数集合，\(p_j\) 为第 \(j\) 个质数。
设 \(g(n,j)\) 表示 \(\sum_{i\leq n}f'(i)[i\in \mathbb{P}\lor\min_{x|i,x\in\mathbb{P}}>p_j]\)。一开始 \(g(n,0)\) 就是 \(\sum_{i\leq n}f'(i)\)，考虑 \(j\) 变大了之后会减少那些数的贡献：
若 \(p_j^2>n\)，则显然 \(g(n,j)=g(n,j-1)\)。
否则，有贡献的是最小质因子恰为 \(p_j\) 的数。可以考虑整体除掉一个 \(p_j\)，因为满足完全积性，所以所有应当删去的树贡献就全部少乘了 \(f'(p_j)\)。
这样就变成了 \(g(\lfloor\frac n {p_j}\rfloor,j-1)\)。但是这样回多算 \([1,\lfloor\frac n {p_j}\rfloor]\) 内质数的贡献，需要再减去。
转移方程为：
\[g(n,j)=\begin{cases} g(n,j-1)\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad p_j^2>n\\ g(n,j-1)-f'(p_j)(g(\lfloor\frac n {p_j}\rfloor,j-1)-\sum_{j\leq \lfloor\frac n {p_j}\rfloor,j\in\mathbb{P}}f'(j))\;\;\;\;\;\;p_j^2\leq n \end{cases} \]这样最后得到的 \(g(n,\lvert\mathbb P\rvert)\) 就是 \(n\) 以内所有质数处的函数值。经过上述操作可以把错误的函数值舍去，最后剩下的是对的函数值。
预处理 \(g\) 之后设 \(S(n,j)\) 表示 \(\sum_{i\leq n}f(i)(\min_{p|i,p\in \mathbb P}p>p_j)\)。
首先对于 \(i\) 是质数的部分，答案就是 \(g(n,\lvert\mathbb P\rvert)\)，这样会多算 \(<j\) 的质数，再减去 \(\sum_{i\in\mathbb P,i<j} f(i)\)。
否则枚举其最小的因子 \(p\) 和次幂 \(e\)，则有：
\[S(n,j)=g(n,\lvert\mathbb P\rvert)-\sum_{i\in\mathbb P,i<j} f(i)+\sum_{i>j}\sum_{e\geq 1}^{p_i^e\leq n}f(p_i^e)(S(\lfloor\frac n {p_i^e}\rfloor,i)+[e>1]) \]加入 \([e>1]\) 是因为质数处的点值已经被计算过，但是 \(p^k(k>1)\) 的点值还没有被计算。
边界条件是 \(S(n,j)=0(p_j>n)\)。但是这样预处理 \(g\) 的复杂度达到了 \(\mathcal O(n\log n)\)。
注意到 \(g(i,j)\) 有用当且仅当存在 \(k\) 使得 \(\lfloor\frac n k\rfloor=i\)，所以有用的点只有 \(\mathcal O(\sqrt n)\) 个。总复杂度是 \(\mathcal O(\mathcal O(\frac{n^{\frac 3 4}}{\log n}))\) 不会证。
ll n,w[300010],g1[300010],g2[300010],sum1[300010],sum2[300010],pr[300010];
int N,m,cnt,id1[300010],id2[300010],iv6;
bool v[300010];
inline void sieve()
{
for(int i=2;i<=N;++i)
{
if(!v[i])pr[++cnt]=i;
for(int j=1;j<=cnt&&pr[j]*i<=N;++j)
{v[pr[j]*i]=1;if(i%pr[j]==0)break;}
}
for(int i=1;i<=cnt;++i)
{
sum1[i]=Cadd(sum1[i-1],pr[i]);
sum2[i]=(sum2[i-1]+1ll*pr[i]*pr[i])%MOD;
}
}
inline int f1(ll x){x%=MOD;return x*(x+1)/2%MOD;}
inline int f2(ll x){x%=MOD;return x*(x+1)%MOD*(2*x+1)%MOD*iv6%MOD;}
inline int getid(ll x){return x<=N?id1[x]:id2[n/x];}
ll S(ll x,int j)
{
if(pr[j]>x)return 0;
ll ans=Cdel(Cdel(g2[getid(x)],g1[getid(x)]),Cdel(sum2[j],sum1[j]));
for(int i=j+1;i<=cnt&&pr[i]*pr[j]<=x;++i)
{
for(ll e=1,sp=pr[i];sp<=x;sp*=pr[i],++e)
ans=(ans+sp%MOD*(sp%MOD-1)%MOD*(S(x/sp,i)+(e>1)))%MOD;
}
return ans;
}
inline void mian()
{
read(n),N=sqrt(n),sieve(),iv6=power(6,MOD-2);
for(ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l),w[++m]=n/l;
g1[m]=f1(w[m])-1,g2[m]=f2(w[m])-1;
if(w[m]<=N)id1[w[m]]=m;else id2[n/w[m]]=m;
}
for(int i=1;i<=cnt;++i)
{
for(int j=1;j<=m&&pr[i]*pr[i]<=w[j];++j)
{
g1[j]=((g1[j]-1ll*pr[i]*(g1[getid(w[j]/pr[i])]-sum1[i-1]))%MOD+MOD)%MOD;
g2[j]=((g2[j]-1ll*pr[i]*pr[i]%MOD*(g2[getid(w[j]/pr[i])]-sum2[i-1]))%MOD+MOD)%MOD;
}
}
cout<<S(n,0)+1;
}