min_25筛学习笔记

min25筛是一个神奇的筛法

是较为常用的筛法

假设要求(f(x))
要求

1. (f)是积性函数
2. (f(p^e))是关于(p)的低阶多项式
3. 存在一个在质数处取值与f一样的完全积性函数

以下以洛谷上的板子题为例
求的积性函数是(f(p^e)=p^e(p^e-1)=(p^e)^2-p^e)

主要思路是分情况讨论,按照实数和合数进行分类

问题变成了求

[sumlimits_{pleq n}f(p)+sumlimits_{i=1}^n[i otin prime]f(i) ]

对后面那个再次进行分类讨论，枚举最小质因子及其指数

因为最小质因子是根号级别的所以可以直接枚举

[sumlimits_{pleq n}f(p)+sumlimits_{pin prime && p^eleq n}f(p^e)(sumlimits_{ileq frac{n}{p^e} && mindiv(i)>p}f(i)) ]

分完情况了我们来进行计算

因为(f(p))是低阶多项式，所以可以拆成单项式，计算(f(p)=p^k)的前缀和

这里就是拆成一次和二次

质数部分

因为

首先肯定是不能进行线性筛的

我们来整一个DP

设(g(n,i)=sumlimits_{j=1}^n[jleq n && mindiv(j)>prime_j]j^k)

设(sp_i=sumlimits_{i=1}^nf(p_i))

用人话说就是所有小于(n)的数中没有小于第(i)个质数的质因子的数

后面的(i^k)是一个和(f)在质数处取值一样的完全积性函数

我们发现(g(n,x))就是最后的答案
其中(p_x)是最后一个小于等于(sqrt{n})的质数

考虑怎样递推

[g(n,i)=g(n,i-1)-p_i^k(g(lfloorfrac{n}{p_i} floor,j-1)-g(p_{i-1},i-1)) ]

[g(n,i)=g(n,i-1)-p_i^k(g(lfloorfrac{n}{p_i} floor,j-1)-sp_{i-1}) ]

考虑这个递推式的含义

从(i-1)到(i)少了一些数
因为是完全积性函数，所以可以直接把最小质因子提出来

后面那个是为了防止把质数给弄没了

因为我们需要的(p)是根号级别的，所以可以直接线性筛暴力算

但是(n)的范围还是太大了，无法都算出来

但我们有结论

[lfloorfrac{lfloorfrac{n}{a} floor}{b} floor=lfloorfrac{n}{ab} floor ]

所以我们的(n)在转移的过程中一直除来除去，也只会用到能表示成(lfloorfrac{n}{x} floor)的值

这显然只有(sqrt{n})个

我们需要在存储dp数组的时候进行一些奥妙重重的实现

发现第二维每次只会增加1，所以可以直接滚掉

计算答案

同样使用dp来计算答案

设(S(n,i))表示最小质因子大于(p_i)的数的函数值之和

那么最终要求得就是(S(n,0))

这里的(g(n))指(g(n,x))，(x)是小于等于(sqrt{n})的最大值

[S(n,i)=g(n)-sp_i+sumlimits_{p_k^eleq n&& k>i)}f(p_k^e)(S(lfloorfrac{n}{p_k^e} floor,k)+[e ot=1]) ]

后面那个中括号的意思差不多是那个数不是质数

然后就做完了

直接爆搜就可以，复杂度是对的

不能处理多组询问

洛谷板子

#include<cstdio>
#include<cmath>
using namespace std;
namespace Patchouli{
  const int N=10000005;
  const int MOD=1000000007;
  const int INV2=500000004;
  const int INV3=333333336;
  long long prime[N],mindiv[N],sp1[N],sp2[N],cnt;
  long long n;
  long long sqr,tot,g1[N],g2[N],w[N],pos1[N],pos2[N];
  inline long long read(){
    long long a=1,b=0;char t;
    do{t=getchar();if(t=='-')a=-1;}while(t>'9'||t<'0');
    do{b=b*10-'0'+t;t=getchar();}while(t>='0'&&t<='9');
    return a*b;
  }
  inline void init(int n){
    mindiv[1]=1;
    for(int i=2;i<=n;++i){
      if(!mindiv[i]){
	mindiv[i]=i;
	prime[++cnt]=i;
	sp1[cnt]=(sp1[cnt-1]+i)%MOD;
	sp2[cnt]=(sp2[cnt-1]+1ll*i*i)%MOD;
      }
      for(int j=1;j<=cnt&&i*prime[j]<=n&&prime[j]<=mindiv[i];++j)
	mindiv[i*prime[j]]=prime[j];
    }
  }
  long long S(long long x,int y){
    if(prime[y]>=x)return 0;
    long long k=x<=sqr?pos1[x]:pos2[n/x];
    long long ans=(g2[k]-g1[k]+MOD-(sp2[y]-sp1[y])+MOD)%MOD;
    for(int i=y+1;i<=cnt&&prime[i]*prime[i]<=x;++i){
      long long tmp=prime[i];
      for(int e=1;tmp<=x;++e,tmp=tmp*prime[i]){
	long long res=tmp%MOD;
	ans=(ans+res*(res-1)%MOD*(S(x/tmp,i)+(e!=1)))%MOD;
      }
    }
    return ans%MOD;
  }
  int QAQ(){
    n=read();
    sqr=sqrt(n);
    init(sqr);
    for(long long l=1,r;l<=n;l=r+1){
      r=n/(n/l);
      w[++tot]=n/l;
      long long tmp=(n/l)%MOD;
      //求出从2到tmp的自然数之和
      g1[tot]=(tmp*(tmp+1)/2)%MOD-1;
      //求出从2到tmp的平方和，式子写这么奇怪只是为了少取模  
      g2[tot]=tmp*(tmp+1)/2%MOD*(2*tmp+1)%MOD*INV3%MOD-1;
      //记录一下各个值出现的位置
      if(n/l<=sqr)pos1[n/l]=tot;
      else pos2[n/(n/l)]=tot;
    }
    for(int i=1;i<=cnt;++i){
      for(int j=1;j<=tot&&prime[i]*prime[i]<=w[j];++j){
	long long k=w[j]/prime[i];
	k=k<=sqr?pos1[k]:pos2[n/k];
	g1[j]-=prime[i]*(g1[k]-sp1[i-1]+MOD)%MOD;
	g2[j]-=prime[i]*prime[i]%MOD*(g2[k]-sp2[i-1]+MOD)%MOD;
	g1[j]%=MOD;
	g2[j]%=MOD;
      }
    }
    printf("%lld
",(S(n,0)+1+MOD)%MOD);
    return false;
  }
}
int main(){
  return Patchouli::QAQ();
}