一種篩法

這篇文章講的是一種篩法，我我的將它稱之爲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 \rfloor$，求出$\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\n",int(ans));
}

拿DIVCNT2測了測速，大概比以前寫的$O(n^{\frac{2}{3}})$還快一倍= = 固然也多是由於以前那份寫的就醜

upd：修復了一個bug