Orz min25 Orz gtr

Min25 篩藉助一個神奇的類似 DP 陣列的手段降低了複雜度

BTW ，這裡說一句題外話，每次看到 DP 的時候，總是不由得想起在 ATABC226F 題題解中看到的一句話，感到這句話是引人深思，並且蒟蒻現在也沒有參透：（相信大家英文很好）

Before getting down to business, let us state the general rule, not limited to this problem: when trying to reduce the time complexity with DP, it is important to focus on “what aspect can be identified.” It is a very common technique to equate multiple states into one to reduce the computational complexity, enabling us to solve more problems.

讓我們進入正題

擴充套件埃氏篩

學習自與擴充套件埃氏篩（min_25篩？）玩耍

其實感覺擴充套件埃氏篩就是 Min25 的簡約版

我們知道埃氏篩的思想是說把每個素數的倍數全部篩掉，其累贅之處在於一個合數會被多個素數反覆篩到

一個優化思路成就了線性篩，也就是保證每一個合數只會被篩一次

另外一個優化思路是在於有些題我們好像沒必要知道具體有哪些素數

例1 loj6235 區間素數個數

給出一個 \(n\) ，求出 \(1\to n\) 的素數個數，\(n\le 10^{11}\)

設 \(f(x)\) 表示 \(1\dots x\) 之間的素數個數

考慮最開始我們只確定了 \(1\) 不是素數，那麼有 \(f(x)=x-1\)

假設現在我們處理到了素數 \(p\) ，即需要把所有包含 \(p\) 的合數篩掉

那麼，對於 \(f(x)\) 而言，\(1\dots x\) 中，有 \(\frac xp\) 個數是 \(p\) 的倍數，考慮到直接 \(f(x)-\frac xp\)​ 的話可能會算重，所以我們每次應該幹掉的是以 \(p\) 為最小質因子的合數（線性篩的思想）。

此時的 \(f(\frac xp)\) 存的是什麼呢？存的是所有的素數和最小質因子 \(\ge p\) 的合數

那我們直接 \(f(x)-=f(\frac xp)\) 就可以了？不行，我們損失了所有\(< p\) 的素數，那麼最後得到了 \(f(x)-=f(\frac xp)-f(p-1)\)

​

埃氏篩只需要處理 \(\sqrt n\) 內的素數

考慮最後答案是 \(f(n)\) ，\(f(n)\) 需要用到的是 \(f(P)\) 和 \(f(\frac nP)\) ，其中 \(P\) 表示所有素數。

顯然， \(\frac np\) 的取值只有 \(\sqrt n\)​​ 種，所以我們記 \(f2(x)\) 表示 \(1\dots\frac nx\) 之間的素數，這樣 \(f,f2\) 只需要開 \(\sqrt n\) 空間，也就是使用 \(f(x)\) 記錄 \(1\dots\sqrt n\) ，\(f2(x)\) 記錄 \(\sqrt n\dots n\)

二者互相推求即可，答案 \(f2(1)\)​ ，據說複雜度 \(O(\dfrac{n^{\frac 34}}{ln\ n})\)

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6+5;

int read(){
    int s=0,w=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')s=s*10+ch-'0',ch=getchar();
    return s*w;
}

int x,f1[N],f2[N];

signed main(){
	x=read(),lim=sqrt(x);
    for(int i=1;i<=lim;++i)f1[i]=sum(i),f2[i]=sum(x/i);
    for(int i=2;i<=lim;++i)if(f1[i]^f1[i-1]){// i 是素數 
        for(int j=1;j<=lim/i;++j)f2[j]-=f2[j*i]-f1[i-1];
        for(int j=lim/i+1;j<=x/i/i&&j<=lim;++j)f2[j]-=f1[x/j/i]-f1[i-1];
        //這個 x/j/i 怎麼來的可以把 j=lim/i+1 帶進去，容易發現和 f2[j] 意義相同 
        for(int j=lim;j>=i*i;--j)f1[j]-=f1[j/i]-f1[i-1];//從大到小更新 f1，防止重複
    }
	printf("%lld\n",f2[1]); 
    return 0;
}

類似題：loj6202 葉氏篩法

Min25 篩

此版塊下的所有 \(p\) 預設為質數，\(p_i\) 表示第 \(i\) 個質數，比如 \(p_3=5\) 。定義 \(lpf(i)\) 為 \(i\) 的最小質因數（\(\text{Lowest Prime Factor}\)）。

假設我們想求 \(S(n)=∑_{i≤n}f(i)\)

Min25 篩適用於一些積性函式 \(f\)，其中 \(f(p)\) 可以用多項式表示，且 \(f(p^k)\) 很好求值。

對於埃氏篩，我們是每次將最小質因子為 \(p\) 的合數篩去，剩下的就是質數。

我們知道合數最小質因子至多為 \(\sqrt n\)，所以合數的貢獻可以通過列舉最小質因子來計算，質數我們則使用另外的方法。

也就是說，我們把 \(S(n)\) 拆成這樣：

\[S(n)=\sum_{i\in P}f(i)+\sum_{i\in P,i\leq\sqrt n,i^e\leq n}f(i^e)\sum_{1\leq x\leq n/i^e,lpf(x)> i} f(x) \]

即是素數的貢獻+以素數 \(i\) 為最小質因子的合數的貢獻

素數貢獻

假裝現在我們正在解決這道例題： P5325 【模板】Min_25篩 - 洛谷

為了計算素數貢獻，Min25 篩先把 \(f\) 看成完全積性函式 \(g\)，其中 \(f(x),g(x)\) 在 \(x\in P\) 處相同。這也就是為什麼要求 \(f(p)\) 可以用多項式表示的原因：多項式可以被拆成單項式，而對單項式而言計算 \(g\) 是容易的。

——GTR

因為 \(f(p^k)=p^k(p^k-1)=p^{2k}-p^k\) ，所以我們不妨平方項和一次向單獨看成兩個多項式處理

令 \(g(p)=p^k,g(x)=\Pi f(p_i^a)\) ,顯然，當 \(x\notin P\) 時，\(g(x)\neq f(x)\)

所謂“單項式而言，計算 \(g\) 是容易的，也就是一點小學數學知識：\(\sum_{i=1}^ni=\frac{n(n+1)}2,\sum_{i=1}^ni^2=\frac{n(n+1)(2n+1)}6,\sum_{i=1}^ni^3=(\frac{n(n+1)}2)^2\)

開始計算：

我們可以從前面擴充套件埃氏篩那裡毛一個“DP”陣列過來：

設 \(G(n,j)\) 表示前 \(n\) 個數，經過前 \(j\) 個質數的篩選，剩下的數的貢獻

標準化的，\(G(n,j)=\sum_{i\leq n}g(i)[i\in P \or lpf(i)>p_j]\)

所謂剩下的數，就是指的是 \(lpf>j\) 的合數和素數，大佬們可以自行理解一下這個 \(G\) 和擴充套件埃氏篩裡面那個 \(f\) 的等價關係

顯然 \(G(n,0)=\sum_{i\leq n}g(i)\)

顯然素數的貢獻就是 $ G(n,|P|)$ ，我們考慮遞推這個 \(g\)，發現 \(G(n,j-1)\) 比 \(G(n,j)\) 多的是

\[G(n,j-1)-G(n,j)=\sum_{i\leq n}g(i)[i\notin P\and lpf(i)=p_j] \]

考慮從小到大遞推，辣麼有：

\[G(n,j-1)-G(n,j)=\begin{cases} 0,&p_j^2>n\\ g(p_j)(G(\frac n{p_j},j-1)-G(p_j,j-1)),&p_j^2\leq n \end{cases} \]

怎麼理解呢，埃氏篩的式子叫做：\(f(x)-=f(\frac xp)-f(p-1)\)

而我們這裡僅僅是多加了一維（反正埃氏篩定義的 \(f\) 其實也是滾動的），而且多算上一個 \(g(p_j)\) 的貢獻

即是，提出一個最小質因子 \(p_j\) 後，那麼在剩下的數中就不能有小於它的質因子了，也就是 \(G(\frac n{p_j},j-1)\) ，但是這樣多減掉了一些素數，哪些？小於 \(p_j\) 的所有素數的貢獻都沒了

所以可以發現 \(G(p_j,j-1)=\sum_{i\leq j}g(p_i)\)

顯然，我們沒有必要求出每一個 \(G(i,x)\) ，有一個結論叫做：\(\lfloor\dfrac{\lfloor\frac na\rfloor}{b}\rfloor=\lfloor\dfrac n{ab}\rfloor\)

所以我們只需要求出 \(\lfloor\frac nx\rfloor\) 形式的數就可以了，這樣的數有 \(\sqrt n\) 個，仿照埃氏篩 \(f2\) 式的處理，我們如此離散化:

可以用 \(ind1[x]\) 表示 \(x\) 這個數對應的陣列下標，\(ind2[x]\)表示 \(n/x\) 這個數對應的下標。這樣兩個 \(ind\) 陣列最大都只會到 \(\sqrt n\)。

複雜度 O(能過) \(O(\dfrac{n^{3/4}}{\log n})\) ，證明可以看 這裡

總貢獻

注意這裡我們把素數和合數貢獻揉在一起算了

rua，不同地，我們記 \(F(n,j)\) 表示 \(\sum_{i\leq n}f(i)[lpf(i)>j]\) ，答案 \(F(n,0)+f(1)\)

\[F(n,j)=G(n,|P|)-G(p_j,j-1)+\sum_{i>j,p_i^e\leq n}f(p_i^e)(F(\frac n{p_i^e},i)+[e\neq 1]) \]

即是大於 \(p_j\) 的質數，然後列舉最小質因子 \(p_i>p_j\) ，\([e=1]\) 時，就是素數 \(p_i\) 貢獻，算過了。

遞迴求解即可，總複雜度 \(O(n^{0,75}/ln\ n+n/poly(ln\ n))\)

例題程式碼

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6+5;
const int mod = 1e9+7;

char buf[1<<23],*p1=buf,*p2=buf;
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
int read(){
    int s=0,w=1; char ch=getchar();
    while(!isdigit(ch)){if(ch=='-')w=-1;ch=getchar();}
    while(isdigit(ch))s=s*10+(ch^48),ch=getchar();
    return s*w;
}

int prime[N],sp1[N],sp2[N];
int w[N],g1[N],g2[N],ind1[N],ind2[N];
int n,sqr,cnt,tot;
bool vis[N];

void getprime(int n){
    for(int i=2;i<=n;++i){
        if(!vis[i])
            prime[++cnt]=i,
            sp1[cnt]=(sp1[cnt-1]+i)%mod,
            sp2[cnt]=(sp2[cnt-1]+i*i%mod)%mod;
        for(int j=1;j<=cnt&&prime[j]*i<=n;++j){
            vis[prime[j]*i]=true;
            if(i%prime[j]==0)break;
        }
    }
}
int calc(int x,int y){
    if(prime[y]>=x)return 0;
    int k=x<=sqr?ind1[x]:ind2[n/x];
    int res=(g2[k]-g1[k]+mod-sp2[y]+sp1[y]+mod)%mod;
    for(int i=y+1;i<=cnt&&prime[i]*prime[i]<=x;++i){
        int mp=prime[i],t=mp%mod;
        for(int j=1;mp<=x;++j,mp=mp*prime[i],t=mp%mod)
            (res+=t*(t-1)%mod*(calc(x/mp,i)+(j!=1))%mod)%=mod;
    }
    return res;
}

signed main(){
    n=read(),sqr=sqrt(n),getprime(sqr);
    for(int i=1;i<=n;){
        int j=n/(n/i);
        w[++tot]=n/i,g1[tot]=w[tot]%mod;
        g2[tot]=(g1[tot]*(g1[tot]+1)%mod*(2*g1[tot]+1)%mod*(166666668)%mod-1+mod)%mod;
        g1[tot]=(g1[tot]*(g1[tot]+1)/2%mod-1+mod)%mod;
        w[tot]<=sqr?ind1[w[tot]]=tot:ind2[n/w[tot]]=tot;
        i=j+1;
    }
    for(int i=1;i<=cnt;++i)
        for(int j=1;j<=tot&&prime[i]*prime[i]<=w[j];++j){
            int k=w[j]/prime[i];
            k=k<=sqr?ind1[k]:ind2[n/k];
            g1[j]-=prime[i]*(g1[k]-sp1[i-1]+mod)%mod;
            g2[j]-=prime[i]*prime[i]%mod*(g2[k]-sp2[i-1]+mod)%mod;
            (g1[j]+=mod)%=mod,(g2[j]+=mod)%=mod;
        }
    printf("%lld\n",(calc(n,0)+1)%mod);
    return 0;
}

Min25 應用

我不會，等我長大後再學習 Min25篩應用_一位蒟蒻的小部落格-CSDN部落格_min25篩的用途