這個演算法可以認為是洲閣篩的一種簡易替代品(另一種寫法),時空常數上和程式碼複雜度上都比洲閣篩優秀,時間複雜度基本和洲閣篩是一樣的.
使用條件和洲閣篩一樣,對積性函式求和.下面是具體流程

假設F(x)為積性函式,且如果x為質數,F(x)=∑xk$F(x)=\sum x^k$
按照洲閣篩的方法,對於每個k$k$分開統計貢獻
對於每個不超過n−−√$\sqrt{n}$的質數p$p$,設g(x)=∑[x是質數或x的每個質因子都不小於p]xk$g(x)=\sum [x是質數或x的每個質因子都不小於p]x^k$
從小到大列舉p$p$,再從大到小列舉所有要求的x>=p2$x>=p^2$,令g(x)−=(g(x/p)−g(p−1))∗pk$g(x)$

再定義S(n,x)代表∑ni=2[i的最小質因子不小於x]F(i)$\sum _{i=2}^n[i的最小質因子不小於x]F(i)$
求S的時候現將質數的貢獻加進來,這個可以通過g算出來
接下來從p開始往上列舉質數p
如果p2$p^2$大於n$n$就break,否則列舉正整數j$j$,滿足pj+1<=n$p^{j+1}<=n$,把S(n/pj,p的下一個質數)∗F(pj)+F(pj+1)$S(n/p^j,p的下一個質數)\ast 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;
}