一种筛法

这篇文章讲的是一种筛法，我个人将它称之为Min_25筛。

它可以用来求积性函数$F(x)$的前缀和，条件与洲阁筛一样，可以快速地对一段质数的F求和。

它可以替代洲阁筛，而且空间常数、时间常数、代码复杂度远比洲阁筛优秀，甚至可以与杜教筛相媲美

时间复杂度与洲阁筛相同据说就是个好写点的洲阁筛

参考链接：

https://post.icpc-camp.org/d/782-spoj-divcnt3/2

https://loj.ac/submission/56015

https://github.com/zimpha/competitive-programming/blob/a0d1ea23778561d29b1d9e4c95eeea63e4e9775a/zoj/3808.cc

首先我们先考虑洲阁筛里面干了什么，首先我们需要对每一个$x=lfloor n/i floor$，求出$sum_{i=1}^x [i是质数] i^k$。

还是类似洲阁筛，我们考虑每一个$leq sqrt{n}$的质数p，对每个x维护$A(x)=sum_{i=1}^x [i是质数或i的每个质因子都>p] i^k$。

我们考虑从大到小更新，我们发现我们要做的事情实际上是对于每个$x geq p^2$，令$A(x)-=(A(x/p)-A(p-1))*p^k$。直接暴力更新。这个部分和洲阁筛复杂度是一样的，也是$O(frac{n^{frac{3}{4}}}{log(n)})$的。

接下来考虑如何对函数求和，我们可以十分暴力地求和。我们定义S(n,x)表示对于2到n的，每个质因子都$geq x$的数的F之和。

考虑如何求S。我们先对质数的F求和，这个用预处理好的A就行。

否则我们从x开始枚举数的最小质因子p，如果$p^2>n$那么break，否则的话我们枚举一个满足$p^{e+1} leq n$的正整数e，把$S(n/p^e,p的下一个质数)F(p^e)+F(p^{e+1})$计入答案。

这部分的复杂度据说也和洲阁筛相同。

UPD：这部分复杂度实际上比较玄学，事实上它的复杂度是 $Theta(n^{1-omega})$ 的，但是对于 $n leq 10^{13}$ 这样的数据范围还是很快的。证明可以参见2018年集训队论文 朱震霆《一些特殊的数论函数求和问题》。

有一些常数优化技巧详见代码。（loj6053）

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <string>
#include <bitset>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <algorithm>
#include <sstream>
#include <stack>
#include <iomanip>
using namespace std;
#define pb push_back
#define mp make_pair
typedef pair<int,int> pii;
typedef long long ll;
typedef double ld;
typedef vector<int> vi;
#define fi first
#define se second
#define fe first
#define FO(x) {freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);}
#define Edg int M=0,fst[SZ],vb[SZ],nxt[SZ];void ad_de(int a,int b){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;}void adde(int a,int b){ad_de(a,b);ad_de(b,a);}
#define Edgc int M=0,fst[SZ],vb[SZ],nxt[SZ],vc[SZ];void ad_de(int a,int b,int c){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;vc[M]=c;}void adde(int a,int b,int c){ad_de(a,b,c);ad_de(b,a,c);}
#define es(x,e) (int e=fst[x];e;e=nxt[e])
#define esb(x,e,b) (int e=fst[x],b=vb[e];e;e=nxt[e],b=vb[e])
typedef unsigned us;
typedef unsigned long long ull;
const us MOD=1e9+7;
#define SS 2333333
ll n,c0[SS],c1[SS],b0[SS],b1[SS]; int S,ps[SS/10],pn=0;
inline ull F(ull x,us g)
{
    if(x<=1||ps[g]>x) return 0;
    ull ans=((x>S)?b1[n/x]:c1[x])-c1[ps[g-1]]+MOD; //注意这里原来有bug
    for(us j=g;j<=pn&&ps[j]*(ll)ps[j]<=x;++j)
    {
        ull cn=x/ps[j],ce=ps[j]*(ll)ps[j];
        for(us e=1;cn>=ps[j];++e,cn/=ps[j],ce*=ps[j])
            ans+=F(cn,j+1)*(ps[j]^e)+(ps[j]^(e+1)),ans%=MOD;
    }
    return ans%MOD;
}
int main()
{
    cin>>n; S=sqrtl(n);
    for(int i=1;i<=S;++i)
    {
        ll t=(n/i)%MOD; b0[i]=t;
        b1[i]=t*(t+1)/2%MOD; c0[i]=i;
        c1[i]=i*(ll)(i+1)/2%MOD;
    }
    for(int i=2;i<=S;++i)
    {
        if(c0[i]==c0[i-1]) continue; //not a prime
        ll x0=c0[i-1],x1=c1[i-1],r=(ll)i*i; ps[++pn]=i;
        int u=min((ll)S,n/(i*(ll)i)),uu=min(u,S/i);
        for(int j=1;j<=uu;++j)
            b1[j]-=(b1[j*i]-x1)*i,
            b0[j]-=b0[j*i]-x0;
        ll t=n/i;
        if(t<=2147483647)
        {
        int tt=t;
        for(int j=uu+1;j<=u;++j)
            b1[j]-=(c1[tt/j]-x1)*i,
            b0[j]-=c0[tt/j]-x0;
        }
        else
        {
        for(int j=uu+1;j<=u;++j)
            b1[j]-=(c1[t/j]-x1)*i,
            b0[j]-=c0[t/j]-x0;
        }
        for(int j=S;j>=r;--j)
            c1[j]-=(c1[j/i]-x1)*i,
            c0[j]-=c0[j/i]-x0;
    }
    for(int i=1;i<=S;++i)
    {
        c1[i]-=c0[i];
        b1[i]-=b0[i];
        if(i>=2) c1[i]+=2;
        if(n>=2LL*i) b1[i]+=2;
        c1[i]=(c1[i]%MOD+MOD)%MOD;
        b1[i]=(b1[i]%MOD+MOD)%MOD;
    }
    ps[pn+1]=S+1;
    ll ans=1+F(n,1);
    ans=(ans%MOD+MOD)%MOD;
    printf("%d
",int(ans));
}

拿DIVCNT2测了测速，大概比之前写的$O(n^{frac{2}{3}})$还快一倍= = 当然也可能是因为之前那份写的就丑

upd：修复了一个bug