洲阁筛

定义非完全积性函数F(x)
我们现在要做的事情是求：

∑x=1nF(x)

还要满足的是：F(pc)是与p相关的低阶多项式
考虑将每个x的质因数分成小于等于n√和大于n√两部分：

∑x=1nF(x)=∑1≤i≤ni没有大于n√的质因数F(i)⎛⎝⎜⎜⎜1+∑n√<j≤⌊ni⌋j是质数F(j)⎞⎠⎟⎟⎟

注意到当i>n√时后面部分的贡献是1，所以要处理的东西有两部分：

∑1≤i≤n√∑n√<j≤⌊ni⌋j是质数F(j)

和

∑n√<i≤ni没有大于n√的质因数F(i)

先解决前面的部分：
设gk(i,j)表示[1,j]内与前i个质数互质的数的k次方和：
那么有转移：

gk(i,j)=gk(i−1,j)−pik×gk(i−1,⌊jpi⌋)

考虑gk(i,j)，当pi+1>j时gk(i,j)的值为1，即当pi>⌊jpi⌋时，后面的贡献为pik
就是说，每次处理gk(i,j)的时候，从小到大枚举i，如果当前的pi2>j，就直接跳出，而调用这个位置之后的g时就直接预处理一下质数的k次方前缀和
由于n以内的质数密度是O(1logn)的，j的取值只有n√种
所以复杂度是这样分析的：

∑n√i=1(ni−−√+i√)logn

分子的部分为：

∫n√0x12dx=23(n12)32=23n34

和

∫n√0nx−−√dx=n12∫n√0x−12dx=n12×2(n12)12=2n34

所以这部分的时间复杂度是O⎛⎝n34logn⎞⎠
接着考虑处理后面的部分，回顾一下：

∑n√<i≤ni没有大于n√的质因数F(i)

设f(i,j)表示[1,j]中仅包含小于等于n√的后i个质数的数的F(x)的和
那么有转移：

f(i,j)=f(i−1,j)+∑c≥1F(pic)⋅f(i−1,⌊jpic⌋)

对于f(i,j)，当pi>j时f(i,j)=F(1)=1，就是说后面的部分当pi2>j时后面的贡献是F(pi)
然后就类似前面计算g，一样做就好了
这部分的复杂度跟前面类似
撒花

代码：
待填坑

#include<iostream>#include<cstring>#include<cstdio>#include<algorithm>#include<cmath>#include<set>#include<bitset>#include<map>#define fo(i,a,b) for(int i=a;i<=b;i++)#define fd(i,a,b) for(int i=a;i>=b;i--)using namespace std;typedef long long LL;typedef double db;const int N = 200010;const int L = 16000000;const int mo = 1e+9+7;const int K = 1;int st[N],id[N],s;LL pri[K+1][N],pre[K+1][N];int k;bool bz[N];LL val[N];LL g[K+1][L],f[L],gs[K+1];int tr[L];LL n,qt;LL a[K+1],F[N];int m;map<LL,int>mp;int lim[N];LL add(LL x,LL y){    return x+y>=mo?x+y-mo:x+y;}LL del(LL x,LL y){    return x>=y?x-y:x+mo-y;}void prepare(){    for(LL i=1;i<=n;i=n/(n/i)+1)val[++m]=n/i;    qt=sqrt(n);    F[1]=1;    fo(i,2,qt){        if (!F[i]){            pri[1][++k]=i;            F[i]=i-1;        }        fo(j,1,k){            if (i*pri[1][j]>qt)break;            if (i%pri[1][j]==0){                F[i*pri[1][j]]=F[i]*pri[1][j]%mo;                break;            }            F[i*pri[1][j]]=F[i]*(pri[1][j]-1)%mo;        }    }    fo(i,1,m/2)swap(val[i],val[m-i+1]);    fo(i,1,k){        LL cnt=1;        fo(j,0,K){            pri[j][i]=cnt;            pre[j][i]=add(pre[j][i-1],cnt);            cnt=cnt*pri[1][i]%mo;        }    }    s=0;    fo(i,1,m){        lim[i]=lim[i-1];        while(lim[i]<k&&pri[1][lim[i]+1]<=val[i])lim[i]++;        mp[val[i]]=i;        st[i]=s+1;        id[i]=k;        fo(j,1,k)        if (1ll*pri[1][j]*pri[1][j]>val[i]){id[i]=j-1;break;}        s+=id[i]+1;        tr[st[i]]=i;        fo(j,1,id[i])tr[st[i]+j]=mp[val[i]/pri[1][j]];    }    //----------------------    a[0]=mo-1;    a[1]=1;}LL getf(int i,int j){    if (j>lim[i])return 1;    if (j>id[i])return (((pre[1][lim[i]]+mo-pre[1][j-1])%mo+mo-(lim[i]-j+1))%mo+1)%mo;    return f[st[i]+j];}LL solve(){    //calc g    fo(i,1,m){        //----------------        g[0][st[i]]=val[i]%mo;        if (val[i]%2)g[1][st[i]]=(val[i]+1)/2%mo*(val[i]%mo)%mo;        else g[1][st[i]]=val[i]/2%mo*(val[i]%mo+1)%mo;        int u=i;        fo(j,1,id[i]){            while(val[u]>val[i]/pri[1][j])u--;            if (j-1<=id[u])                fo(x,0,K)                g[x][st[i]+j]=del(g[x][st[i]+j-1],pri[x][j]*g[x][st[u]+j-1]);            else                fo(x,0,K)                g[x][st[i]+j]=del(g[x][st[i]+j-1],pri[x][j]*(del(g[x][st[u]+id[u]],del(pre[x][j-1],pre[x][id[u]]))));        }    }    int u=m;    LL ans=0;    fo(i,1,qt){        while(val[u]>n/i)u--;        fo(x,0,K){            LL tmp=del(g[x][st[u]+id[u]],add(del(pre[x][k],pre[x][id[u]]),1));            ans=(ans+(tmp*a[x]%mo)*F[i]%mo)%mo;        }        //ans=add(ans,F[i]);    }    fo(i,1,m){        LL cnt=del(del(pre[1][lim[i]],pre[1][id[i]]),(lim[i]-id[i]))+1;        fd(j,id[i],1){            LL x=val[i];            LL tmp=pri[1][j]-1;            if (j==id[i])f[st[i]+j]=cnt;            else f[st[i]+j]=f[st[i]+j+1];            for(int now=i;x>=pri[1][j];){                if (x>=pri[1][j]*pri[1][j]){                    now=tr[st[now]+j];                    f[st[i]+j]=add(f[st[i]+j],tmp*getf(now,j+1)%mo);                }                else f[st[i]+j]=add(f[st[i]+j],tmp);                x=x/pri[1][j];                tmp=tmp*pri[1][j]%mo;            }        }    }    ans=add(ans,f[st[m]+1]);    return ans;}int main(){    freopen("data.in","r",stdin);    scanf("%lld",&n);    prepare();    LL ans=solve();    printf("%lld\n",ans);    return 0;}

不知道为什么跑不过10^10….要找张贡提升一下知识水平啦2333