洲閣篩詳解

你還真信了

丟鏈接

這篩對積性函數的要求不同於杜教篩，只消函數在自變量為質數或質數整數冪時是一個低階多項式即可。以下n<=1e11。

首先有一個性質：1~n的每個數，大於$\sqrt{n}$的質因子只有一個。根據是否有大於$\sqrt{n}$的質因子，再根據他是積性函數，得

$\sum_{i=1}^nf(i)=\sum_{i=1}^{n}[i \ \ has \ \ no \ \ prime \ \ greater \ \ than \ \ \sqrt{n}]f(i)(1+\sum_{i=\sqrt{n}}^{\left \lfloor \frac{n}{i} \right \rfloor}[i \ \ is \ \ prime]f(i))$

好的！來觀察他。首先$\sum_{i=1}^{\sqrt{n}}[i \ \ has \ \ no \ \ prime \ \ greater \ \ than \ \ \sqrt{n}]$是可以線性篩的。

於是我們就差兩部分：Part1：$\sum_{i=\sqrt{n}}^{\left \lfloor \frac{n}{i} \right \rfloor}[i \ \ is \ \ prime]f(i)$，Part2：$\sum_{i=\sqrt{n}}^n[i \ \ has \ \ no \ \ prime \ \ greater \ \ than \ \ \sqrt{n}]f(i)$，並希望用$\leq \sqrt{n}$的質數來處理他們。

Part1：$g(i,j)$--$[1,j]$中與前$i$個質數互質的數的k次冪和。有遞推式$g(i,j)=g(i-1,j)-p_i^kg(i-1,\left \lfloor \frac{j}{p_i} \right \rfloor)$。直接算的復雜度是$\frac{n}{log(n)}$。

注意到$p_{i+1}>j$時$g(i,j)=1$，緊接着$p_i>\left \lfloor \frac{j}{p_i} \right \rfloor$即$p_i^2>j$時$g(i-1,\left \lfloor \frac{j}{p_i} \right \rfloor)=1$，因此$g(i,j)=g(i-1,j)-p_i^k$。

也就是說算到$p_i^2>j$時就不必再算了，后面要用某個$g(i,j)$時，若$p_{i+1}>j$則$g(i,j)=1$，否則用$g(i0_j,j)$減去一段$p_i^k$即可，其中$i0_j$表示最慢枚舉到$j$的$i$。

Part2：$h(i,j)$--$[1,j]$中用前$i$個質數湊起來的數的$f$值的和。有$h(i,j)=h(i-1,j)+\sum_{c=1}^{p_i^c<=j}f(p_i^c)h(i-1,\left \lfloor \frac{j}{p_i^c} \right \rfloor)$

然而這里沒有上面那個用平方來卡$j$上限的性質！！咋整！！

反過來。把前$i$個改成后$i$個，注意是$\leq \sqrt{n}$的后$i$個。遞推式同理。但！

$p_i>j$時$h(i,j)=1$，$p_i^2>j$時用后$i-1$個質數湊得$\leq j$的部分只有1，因此$h(i,j)=h(i-1,j)+f(p_i)$。棒！用類似方法解決。

注意$h$實際上回答了$\sum_{i=1}^{n}[i \ \ has \ \ no \ \ prime \ \ greater \ \ than \ \ \sqrt{n}]f(i)$，因此最后加上即可，括號里那個1就不必理了。然后$g$中都是包含1的，但是我不要，因此在算答案的時候記得-1。

以下是那道萬年不變的例題的代碼。

//#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
//#include<map>
#include<math.h>
//#include<time.h>
//#include<complex>
#include<algorithm>
using namespace std;

#define LL long long
int T; LL n;
int m=320000;
#define maxn 640011
LL prime[maxn],f[maxn],sf[maxn],s[maxn],i0[maxn],who[maxn],g[maxn],ff[maxn],mi[maxn],ci[maxn]; 
int lp,lw; bool notprime[maxn];

#define maxh 1000007
struct Hash
{
    int first[maxh],le; struct Edge{LL to,v; int next;}edge[maxn];
    void clear(LL n) {le=2; for (LL i=1;i<=n;i=n/(n/i)+1) first[n/i%maxh]=0;}
    void insert(LL x,LL y) {int z=x%maxh; Edge &e=edge[le]; e.to=x; e.v=y; e.next=first[z]; first[z]=le++;}
    LL find(LL x) {int z=x%maxh; for (int i=first[z];i;i=edge[i].next) if (edge[i].to==x) return edge[i].v; return -1;}
}h;

void pre(int &n)
{
    lp=0; f[1]=sf[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!notprime[i]) {prime[++lp]=i; f[i]=4; mi[i]=i; ci[i]=1;}
        sf[i]=sf[i-1]+f[i]; s[i]=s[i-1]+!notprime[i];
        for (int j=1,tmp;j<=lp && 1ll*i*prime[j]<=n;j++)
        {
            notprime[tmp=i*prime[j]]=1;
            if (i%prime[j]) {f[tmp]=f[i]*4; mi[tmp]=prime[j]; ci[tmp]=1;}
            else {f[tmp]=f[i/mi[i]]*(3*ci[i]+4); mi[tmp]=mi[i]*prime[j]; ci[tmp]=ci[i]+1; break;}
        }
    }
    lp--; n=prime[lp+1]-1;
}

void mg()
{
    lw=0; h.clear(n);
    for (LL i=1;i<=n;i=n/(n/i)+1) lw++,who[lw]=g[lw]=n/i,h.insert(who[lw],lw);
    memset(i0,0,sizeof(i0));
    for (int i=1;i<=lp;i++)
        for (int j=1;j<=lw && prime[i]*prime[i]<=who[j];j++)
        {
            int k=h.find(who[j]/prime[i]);
            g[j]-=g[k]-(i-1-i0[k]);
            i0[j]=i;
        }
}

void mff()
{
    for (int j=1;j<=lw;j++) ff[j]=1;
    for (int i=lp;i;i--)
        for (int j=1;j<=lw && 1ll*prime[i]*prime[i]<=who[j];j++)
        {
            if (prime[i+1]>who[j]) ff[j]=1;
            else if (prime[i+1]*prime[i+1]>who[j]) ff[j]=(s[min((LL)m,who[j])]-s[prime[i+1]-1])*4+1;
            LL tmp=prime[i],k2;
            for (int c=1;tmp<=who[j];tmp*=prime[i],c++)
            {
                LL k=who[j]/tmp;
                if (prime[i+1]>k) k2=1;
                else if (prime[i+1]*prime[i+1]>k) k2=(s[min((LL)m,k)]-s[prime[i+1]-1])*4+1;
                else k2=ff[h.find(k)];
                ff[j]+=k2*(3*c+1);
            }
        }
}

int main()
{
    scanf("%d",&T);
    pre(m);
    while (T--)
    {
        scanf("%lld",&n);
        if (n<=m) {printf("%lld\n",sf[n]); continue;}
        LL ans=0,last; mg(); mff();
        for (int i=1,id;i<=m;i=last+1)
        {
            LL tmp=n/i,kk;
            if (prime[lp+1]>tmp) kk=1;
            else id=h.find(tmp),kk=g[id]-(lp-i0[id]);
            ans+=(sf[last=min((LL)m,n/tmp)]-sf[i-1])*(kk-1)*4;
        }
        printf("%lld\n",ans+ff[1]);
    }
    return 0;
}

配倆圖，表示$g$和$h$的遞推順序。