Min_25筛学习

从http://www.cnblogs.com/zzqsblog/p/8302815.html习得的一种算法

这个算法可以认为是洲阁筛的一种简易替代品(另一种写法),时空常数上和代码复杂度上都比洲阁筛优秀,时间复杂度基本和洲阁筛是一样的.
使用条件和洲阁筛一样,对积性函数求和.下面是具体流程

假设F(x)为积性函数,且如果x为质数,$F(x)=\sum x^k$
按照洲阁筛的方法,对于每个$k$分开统计贡献
对于每个不超过$\sqrt n$的质数$p$,设$g(x)=\sum [x是质数或x的每个质因子都不小于p]x^k$
从小到大枚举$p$,再从大到小枚举所有要求的$x>=p^2$,令$g(x)-=(g(x/p)-g(p-1))*p^k$
上面这个式子的意义大概是枚举合数中的最小质因子,再把它的贡献减去.这个部分复杂度和洲阁筛是一样的.

再定义S(n,x)代表$\sum_{i=2} ^n[i的最小质因子不小于x]F(i)$
求S的时候现将质数的贡献加进来,这个可以通过g算出来
接下来从p开始往上枚举质数p
如果$p^2$大于$n$就break,否则枚举正整数$j$,满足$p^{j+1}<=n$,把$S(n/p^j,p的下一个质数)*F(p^j)+F(p^{j+1})$算入贡献
这部分的复杂度感觉不好证,不过据说是和洲阁筛差不多的

以SPOJDIVCNT3为例,具体实现见代码.这个做法还是比较容易写错的,而且姿势不好的话常数会大很多.
运行时间上最朴素写法12s,加上优化后4.7s,网上随便找到一个洲阁筛21.9s

#include<bits/stdc++.h>
using namespace std;
#define maxn 2000100
long long Min(long long x,long long y){
    return x<y?x:y;
}
int mysqrt(long long n){
    int x=sqrt(n);
    while(x*(long long)x<n) x++;
    while(x*(long long)x>n) x--;//必须满足S是最大的满足S*S<=n的数
    return x;
}
int S;
long long g[maxn];
int prim[maxn],ilem;
long long n;
void getg(){
    S=mysqrt(n),ilem=0;
    long long *l=g+S,*s=g;
    for(int i=1;i<=S;i++)
        s[i]=i,l[i]=n/i;//s[i]->g(i),l[i]->g(n/i),此时g(i)代表不超过i的质数的数量(包括1)
    for(int p=2;p<=S;p++){
        if(s[p]==s[p-1]) continue;//p不是质数
        long long r=p*(long long)p,v=s[p-1];
        int ed1=S/p,ed2=Min(n/r,(long long)S);
        for(int i=1;i<=ed1;i++)
            l[i]-=l[i*p]-v;
        if(n/p<INT_MAX){
            int m=n/p;
            for(int i=ed1+1;i<=ed2;i++)
                l[i]-=s[m/i]-v;
        }
        else{
            long long m=n/p;
            for(int i=ed1+1;i<=ed2;i++)
                l[i]-=s[m/i]-v;
        }
/*** 如果不这么写的话,常数大概会增加一倍左右.把int和long long分开在n=1e11的数据下有0.4s的优化(开O2) ***/
/*** 论减少除法数量的重要性... ***/
        for(int i=S;i>=r;i--)
            s[i]-=s[i/p]-v;
        prim[++ilem]=p;
    }
    prim[++ilem]=S+1;//这句很容易被遗漏
    for(int i=1;i<=S;i++)
        s[i]=s[i]*4,l[i]=l[i]*4;
}
long long F(long long x,int y){
    if(prim[y]>x) return 0;
    long long res=g[x<=S?x:S+n/x]-g[prim[y]-1],st;
    for(int i=y;i<=ilem;i++){
        st=prim[i];
        if(st*(long long)prim[i]>x) break;
        for(int j=1;;j++){
            if(st*prim[i]>x) break;
            res+=F(x/st,i+1)*(j*3+1)+(j*3+4);
            st=st*prim[i];
        }
    }
    return res;
}
int main(){
//  freopen("DIVCNT3.in","r",stdin);
    int T;
    scanf("%d",&T);
    while(T--){
        scanf("%lld",&n);
        getg();
        printf("%lld\n",F(n,1)+1);//注意把1的贡献算进来.
    }
    return 0;
}