Min-25 篩學習筆記

$\text{By DaiRuiChen007}$

一、簡要介紹

Min-25 篩，是一種能在亞線性時間內求出特定的一類積性函數 $f(i)$ 的字首和的演演算法。

具體來說，Min-25 篩可以在 $\mathcal O(\sqrt n)$ 的空間複雜度與 $\mathcal O(\dfrac{n^{3/4}}{\log n})$ 的時間複雜度內求出 $\sum _{i=1}^n f(i)$ 的值。

一般來說，積性函數 $f(i)$ 要滿足如下幾個條件：

$f(i)$ 在質數處可以被表示成一個較為簡單的多項式。
$f(i)$ 在質數的若干次冪（$p^c$ 處）的點值可以 $\mathcal O(1)$ 計算。

一般來說，在 $n$ 較小（$\le 10^{10}$）時，$\mathcal O(\dfrac{n^{3/4}}{\log n})$ 的複雜度是優於杜教篩的 $\mathcal O(n^{2/3})$ 的。

二、演演算法過程

以下內容均用 $p_i$ 表示第 $i$ 小的質數。

考慮拆分出每個 $n$ 的最小質因子然後暴力計算，但是這樣的質因子是 $\mathcal O(n)$ 級別的，考慮分成合數和質數兩種情況討論，其中合數的最小質因子就會被降到 $\mathcal O(\sqrt n)$ 級別。

由於要根據最小質因子討論，因此記 $S(n,k)$ 表示 $1\sim n$ 中所有最小質因子 $>p_k$ 的 $i$ 的 $f(i)$ 之和。

先不管質數的情況，設 $\sum\limits_{i=1}^{p_i\le n} f(p_i)$ 的值為 $Q(n)$，那麼列舉 $>p_k$ 的最小質因子以及其指數可以得到：

\[S(n,k)=Q(n)-\sum_{i=1}^{i\le k} f(p_i)+\sum_{i=k+1}^{p_i^2\le n}\sum_{c=1}^{p_i^c\le n} f(p_i^c)\times (S(\lfloor n/p_i^c\rfloor,i)+[c>1]) \]

其中 $[c>1]$ 是求 $f(p_i^c)$ 的貢獻，容易發現答案就是 $S(n,0)+f(1)$。

假設 $Q(n)$ 已知，預處理出 $f(p_i)$ 的字首和（容易遞迴過程中發現 $p_k^2\le n$，因此只要處理 $\le \sqrt n$ 的質數），這個式子直接暴力遞迴計算的複雜度就是 $\mathcal O(\dfrac{n^{3/4}}{\log n})$。

接下來考慮計算 $Q(n)$，注意到在遞迴過程中的 $n$ 始終可以表示成某個 $\lfloor n/t\rfloor$ 的值，因此我們只要求出 $Q(n)$ 在這 $2\sqrt n$ 個位置的點值即可。

根據要求，$f(p_i)$ 一定是一個簡單的多項式，我們分次數處理，分別求出 $\sum p_i^0,\sum p_i^1,\sum p_i^2,\dots $ 再整合。以 $\sum p_i^1$ 為例，構造完全積性函數 $f'(i)=i$，容易發現這個函數在質數處與 $f(i)$ 相等，因此我們只要求出 $f'(i)$ 在質數處的和。

設 $G(n,k)$ 表示 $1\sim n$ 中，滿足 $i$ 的最小質因子 $>p_k$ 或 $i$ 是質數的 $f'(i)$ 的和。

考慮 $G(n,k-1)-G(n,k)$ 的值，當 $p_k^2>n$ 時，顯然沒有最小質因子 $=p_k$ 的合數，差為 $0$。

否則，計算的數一定有最小質因子 $p_k$，提出 $f'(p_k)$ 得到：

\[G(n,k-1)-G(n,k)=f'(p_k)\times\left(G(\lfloor n/p_k\rfloor,k-1)-\sum_{i=1}^{k-1}f'(p_i) \right) \]

減去的項是去掉有 $<p_k$ 的質因子的項，因此我們得到 $G(n,k)$ 的遞推式：

\[G(n,k)= \begin{cases} G(n,k-1)&p_k^2>n\\ G(n,k-1)-f'(p_k)\times\left(G(\lfloor n/p_k\rfloor,k-1)-\sum_{i=1}^{k-1}f'(p_i) \right) &p_k^2\le n \end{cases} \]

同樣在遞迴的過程中也只會存取到所有的 $\lfloor n/t\rfloor$ 下標處，且最後 $f'(p_i)$ 的字首和也只要處理到 $\le \sqrt n$ 的所有質數。

邊界條件是 $G(n,0)=n(n+1)/2-1$，$Q(n)=G(n,\infty)$，從小到大列舉 $k$，從大到小列舉 $n$ 轉移即可，複雜度也可以證明為 $\mathcal O(\dfrac{n^{3/4}}{\log n})$。

拓展 - 非遞迴版本的實現

對於 $f(i)$ 的字首和 $S(n)$，我們已經介紹了在 $\mathcal O(\dfrac{n^{3/4}}{\log n})$ 的時間複雜度以及 $\mathcal O(\sqrt n)$ 的空間複雜度內計算 $S(n)$ 的演演算法。

但假如我們想和杜教篩一樣，在複雜度不變的情況下求出 $S(n)$ 在所有 $S(\lfloor n/t\rfloor )$ 處的值，應該怎麼處理呢？

最簡單的想法是在遞迴求 $S(n,k)$ 的時候加上記憶化，這樣確實能保證時間複雜度不退化，但是空間複雜度會變成 $\mathcal O(\dfrac{n^{3/4}}{\log n})$，如果想要空間複雜度是 $\mathcal O(\sqrt n)$ 的做法，需要轉變一下 $S$ 的定義再處理。

注意到在處理 $S(n,k)$ 時最棘手的是形如 $S(n,k)\gets S(n',i)$ 的轉移時 $k$ 與 $i$ 不連續，因此導致沒法簡單遞推求值。

我們認為 $G$ 的遞推式很優秀也是因為轉移式永遠只有 $i=k-1$ 的狀態向 $k$ 轉移，因此考慮讓 $S(n,k)$ 的定義也向 $G(n,k)$ 的定義靠攏，不妨設 $S(n,k)$ 表示 $1\sim n$ 中，滿足 $i$ 的最小質因子 $>p_k$ 或 $i$ 是質數的 $f(i)$ 的和。

依然考慮 $S(n,k-1)-S(n,k)$，在 $p_k^2> n$ 時依然為 $0$，在 $p_k^2\le n$ 時列舉 $p_k$ 的次數得到：

\[S(n,k-1)-S(n,k)=\sum_{c=1}^{p_k^c\le n} f(p_k^c)\times\left(S(\lfloor n/p_k^c\rfloor,k)+[c>1]-\sum_{i=1}^{p_i\le \min(p_k,\lfloor n/p_k^c\rfloor)}f(p_i)\right) \]

但是注意到最後一項求和式並不是很好看，考慮化簡，注意到當 $p_k> \lfloor n/p_k^c\rfloor$ 時，$p_k^{c+1}> n$，那麼 $S(\lfloor n/p_k^c\rfloor,k)$ 只會統計 $\le \lfloor n/p_k^c\rfloor$ 的質數，與求和式正好抵消，因此可以把上式寫成：

\[S(n,k-1)-S(n,k)=\sum_{c=1}^{p_k^{c+1}\le n}f(p_k^{c+1})+ f(p_k^c)\times\left(S(\lfloor n/p_k^c\rfloor,k)-\sum_{i=1}^{k}f(p_i)\right) \]

因此得到 $S(n,k)$ 的遞推式：

\[S(n,k-1)= \begin{cases} S(n,k)&p_k^2>n\\[2ex] S(n,k)+\sum\limits_{c=1}^{p_k^{c+1}\le n}f(p_k^{c+1})+ f(p_k^c)\times\left(S(\lfloor n/p_k^c\rfloor,k)-\sum_{i=1}^{k}f(p_i)\right)&p_k^2\le n \end{cases} \]

從大到小列舉 $p_k$ 再從大到小列舉 $n$ 轉移即可。

其中邊界條件為 $S(n,\infty)=Q(n)$，最終的答案是 $S(\lfloor n/t\rfloor,0)$，時間複雜度 $\mathcal O(\dfrac{n^{3/4}}{\log n})$，空間複雜度 $\mathcal O(\sqrt n)$。

有了這個 Min-25 篩的實現後，我們可以解決形如 $\sum_{i=1}^n f(i)\times g(\lfloor n/i\rfloor)$ 的和式計算問題，其中 $f(i)$ 是積性函數，$g(i)$ 是任意一個關於 $i$ 的函數。

那麼我們根據整除分塊列舉 $\lfloor n/i\rfloor $，則 $g$ 函數的值固定，要求的就是某一段特定區間內 $f(i)$ 的和，容易發現這些區間的端點都是某些 $\lfloor n/t\rfloor$，用上面所介紹的非遞迴型的 Min-25 篩處理即可。

三、程式碼實現

以線性篩模板題為例（Link）下面的程式碼是遞迴的實現：

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=2e5+1,MOD=1e9+7;
bool isco[MAXN];
int p[MAXN],tot;
int val[MAXN]; //需要處理的 n (降序排列)
int idx1[MAXN],idx2[MAXN]; //用於儲存每個可能的 n/t 對應的下標
int g1[MAXN],g2[MAXN]; //p,p^2 在所有 n/t 處的字首和
int fp1[MAXN],fp2[MAXN]; //p,p^2 在所有 <=sqrt(n) 的質數處的字首和
int n,m,B;
inline int idx(int v) {
    //v 對應儲存在陣列裡的下標
    return (v<=B)?idx1[v]:idx2[n/v];
}
inline int S(int n,int k) {
    if(p[k]>n) return 0;
    int ans=((g2[idx(n)]+MOD-g1[idx(n)])%MOD+MOD-(fp2[k]+MOD-fp1[k])%MOD)%MOD;
    for(int i=k+1;i<=tot&&p[i]*p[i]<=n;++i) {
        for(int c=1,v=p[i];v<=n;++c,v*=p[i]) {
            //直接遞迴求值
            ans=(ans+(v-1)%MOD*(v%MOD)%MOD*((c>1)+S(n/v,i))%MOD)%MOD;
        }
    }
    return ans;
}
signed main() {
    scanf("%lld",&n),B=sqrt(n);
    for(int i=2;i<=B;++i) {
        if(!isco[i]) p[++tot]=i;
        for(int j=1;j<=tot&&i*p[j]<=B;++j) {
            isco[i*p[j]]=true;
            if(i%p[j]==0) break;
        }
        //線性篩預處理質數
    }
    for(int i=1;i<=tot;++i) {
        //在 <=sqrt(n) 的質數處的字首和
        fp1[i]=(fp1[i-1]+p[i])%MOD;
        fp2[i]=(fp2[i-1]+p[i]*p[i]%MOD)%MOD;
    }
    for(int l=1,r;l<=n;l=r+1) {
        //整除分塊預處理出點值
        r=n/(n/l),val[++m]=n/l;
        if(val[m]<=B) idx1[val[m]]=m;
        else idx2[n/val[m]]=m;
    }
    for(int i=1;i<=m;++i) {
        int k=val[i]%MOD;
        g1[i]=(MOD+1)/2*k%MOD*(k+1)%MOD;
        g2[i]=(MOD+1)/6*k%MOD*(2*k+1)%MOD*(k+1)%MOD;
        g1[i]=g1[i]?g1[i]-1:MOD-1,g2[i]=g2[i]?g2[i]-1:MOD-1;
        //G(n,0) 的值
    }
    for(int k=1;k<=tot;++k) {
        for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {
            //暴力列舉, 按遞推式計算
            g1[i]=(g1[i]+MOD-p[k]*(g1[idx(val[i]/p[k])]+MOD-fp1[k-1])%MOD)%MOD;
            g2[i]=(g2[i]+MOD-p[k]*p[k]%MOD*(g2[idx(val[i]/p[k])]+MOD-fp2[k-1])%MOD)%MOD;
        }
    }
    printf("%lld\n",(S(n,0)+1)%MOD);
    return 0;
}

下面的程式碼是非遞迴的實現：

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=2e5+1,MOD=1e9+7;
bool isco[MAXN];
int p[MAXN],tot;
int val[MAXN]; 
int idx1[MAXN],idx2[MAXN];
int g1[MAXN],g2[MAXN];
int fp1[MAXN],fp2[MAXN];
int S[MAXN];
int n,m,B;
inline int idx(int v) {
    return (v<=B)?idx1[v]:idx2[n/v];
}
signed main() {
    scanf("%lld",&n),B=sqrt(n);
    for(int i=2;i<=B;++i) {
        if(!isco[i]) p[++tot]=i;
        for(int j=1;j<=tot&&i*p[j]<=B;++j) {
            isco[i*p[j]]=true;
            if(i%p[j]==0) break;
        }
    }
    for(int i=1;i<=tot;++i) {
        fp1[i]=(fp1[i-1]+p[i])%MOD;
        fp2[i]=(fp2[i-1]+p[i]*p[i]%MOD)%MOD;
    }
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l),val[++m]=n/l;
        if(val[m]<=B) idx1[val[m]]=m;
        else idx2[n/val[m]]=m;
    }
    for(int i=1;i<=m;++i) {
        int k=val[i]%MOD;
        g1[i]=(MOD+1)/2*k%MOD*(k+1)%MOD;
        g2[i]=(MOD+1)/6*k%MOD*(2*k+1)%MOD*(k+1)%MOD;
        g1[i]=g1[i]?g1[i]-1:MOD-1,g2[i]=g2[i]?g2[i]-1:MOD-1;
    }
    for(int k=1;k<=tot;++k) {
        for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {
            g1[i]=(g1[i]+MOD-p[k]*(g1[idx(val[i]/p[k])]+MOD-fp1[k-1])%MOD)%MOD;
            g2[i]=(g2[i]+MOD-p[k]*p[k]%MOD*(g2[idx(val[i]/p[k])]+MOD-fp2[k-1])%MOD)%MOD;
        }
    }
    for(int i=1;i<=m;++i) S[i]=(g2[i]+MOD-g1[i])%MOD;
    for(int k=tot;k>=1;--k) {
        for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {
            for(int c=1,v=p[k];v*p[k]<=val[i];) { //根據遞推式計算
                S[i]=(S[i]+(v-1)%MOD*(v%MOD)%MOD*(S[idx(val[i]/v)]+MOD-(fp2[k]-fp1[k]))%MOD)%MOD;
                ++c,v*=p[k],S[i]=(S[i]+(v-1)%MOD*(v%MOD)%MOD)%MOD;
            }
        }
    }
    printf("%lld\n",(S[idx(n)]+1)%MOD);
    return 0;
}

實際執行上，非遞迴版本的實現會慢於遞迴版的實現。

四、複雜度分析

空間複雜度 $\mathcal O(\sqrt n)$，可以參照參考程式碼的實現。

時間複雜度可以看做合法的 $S(n,k) / G(n,k)$ 狀態數，對於每個 $n=\lfloor n/t\rfloor$，合法的 $k$ 共有 $\pi(\sqrt{\lfloor n/t\rfloor})$ 個。

因此複雜度可以估計為：

\[\begin{aligned} \mathcal T(n) &=\sum_{i^2\le n} \mathcal O(\pi(\sqrt i))+\mathcal O(\pi(\sqrt{n/i}))\\[2.5ex] &=\sum_{i^2\le n} \mathcal O\left(\dfrac{\sqrt i}{\ln \sqrt i}\right)+\mathcal O\left(\dfrac{\sqrt{n/i}}{\ln \sqrt{n/i}}\right)\\[2.5ex] &=\mathcal O\left(\dfrac{1}{\log n}\int_1^{\sqrt n} \sqrt{\dfrac nx} \mathrm dx\right)\\[2.5ex] &=\mathcal O\left(\dfrac{\sqrt n}{\log n}\times\left. 2\sqrt x \right|^{\sqrt n}_1 \right)\\[2.5ex] &=\mathcal O(\dfrac{n^{3/4}}{\log n}) \end{aligned} \]

五、例題

I. [SPOJ-34096] DIVCNTK

Problem Link

題目大意

給定 $n$，求 $\sum_{i=1}^n \sigma_0(i^k)$（即 $i^k$ 的因子個數）。

資料範圍：$n\le 10^{10}$。

思路分析

容易證明 $\sigma_0(i)$ 是積性函數，同理 $\sigma_0(i^k)$ 同樣也是積性函數，考慮求 $\sigma_0(p^k)$，顯然 $\sigma_0(p^k)=k+1$，我們只需要求所有 $1\sim \lfloor n/t\rfloor$ 之間的質數個數，構造 $f'(p)=1$，直接用 Min-25 篩處理即可。

時間複雜度 $\mathcal O(\dfrac{n^{3/4}}{\log n})$

程式碼呈現

#include<bits/stdc++.h>
#define int unsigned long long
using namespace std;
const int MAXN=2e5+1;
bool isco[MAXN];
int p[MAXN],tot;
int idx1[MAXN],idx2[MAXN],val[MAXN];
int g[MAXN];
inline void solve() {
    int n,K,m=0,B;
    scanf("%llu%llu",&n,&K),B=sqrt(n);
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l),val[++m]=n/l;
        if(val[m]<=B) idx1[val[m]]=m;
        else idx2[n/val[m]]=m;
        g[m]=val[m]-1;
    }
    auto idx=[&](int v) { return v<=B?idx1[v]:idx2[n/v]; };
    for(int k=1;k<=tot;++k) {
        for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {
            g[i]-=g[idx(val[i]/p[k])]-(k-1);
        }
    }
    auto S=[&](auto self,int n,int k) -> int {
        if(n<=p[k]) return 0;
        int ans=(g[idx(n)]-k)*(K+1);
        for(int i=k+1;i<=tot&&p[i]*p[i]<=n;++i) {
            for(int c=1,v=p[i];v<=n;++c,v*=p[i]) {
                ans+=(c*K+1)*(self(self,n/v,i)+(c>1));
            }
        }
        return ans;
    };
    printf("%llu\n",S(S,n,0)+1);
}
signed main() {
    for(int i=2;i<MAXN;++i) {
        if(!isco[i]) p[++tot]=i;
        for(int j=1;j<=tot&&i*p[j]<MAXN;++j) {
            isco[i*p[j]]=true;
            if(i%p[j]==0) break;
        }
    }
    int T;
    scanf("%llu",&T);
    while(T--) solve();
    return 0;
}

II. [洛谷-4213] SUM

Problem Link

題目大意

求 $\mu(i),\varphi(i)$ 的字首和。

資料範圍：$n\le 2^{31}$。

思路分析

注意到 $\mu(p^c),\varphi(p^c)$ 的計算是容易的，在質數處的點值滿足：$\mu(p)=-1,\varphi(p)=p-1$，因此預處理 $G(n,k)$ 時分別處理 $\sum p_i^0,\sum p_i^1$ 的值即可。

時間複雜度 $\mathcal O(\dfrac{n^{3/4}}{\log n})$。

程式碼呈現

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=2e5+1;
bool isco[MAXN];
int p[MAXN],tot,val[MAXN],fp1[MAXN];
int idx1[MAXN],idx2[MAXN],g0[MAXN],g1[MAXN];
inline void solve() {
    int n,m=0;
    scanf("%lld",&n);
    int B=sqrt(n); tot=0;
    for(int i=2;i<=B;++i) {
        if(!isco[i]) p[++tot]=i,fp1[tot]=i+fp1[tot-1];
        for(int j=1;j<=tot&&p[j]*i<=B;++j) {
            isco[p[j]*i]=true;
            if(i%p[j]==0) break;
        }
    }
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l),val[++m]=n/l;
        if(val[m]<=B) idx1[val[m]]=m;
        else idx2[n/val[m]]=m;
        g0[m]=val[m]-1,g1[m]=val[m]*(val[m]+1)/2-1;
    }
    auto idx=[&](int v) { return (v<=B)?idx1[v]:idx2[n/v]; };
    for(int k=1;k<=tot;++k) {
        for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {
            g0[i]-=g0[idx(val[i]/p[k])]-(k-1);
            g1[i]-=p[k]*(g1[idx(val[i]/p[k])]-fp1[k-1]);
        }
    }
    auto mu=[&](auto self,int n,int k) -> int {
        if(n<=p[k]) return 0;
        int ans=k-g0[idx(n)];
        for(int i=k+1;i<=tot&&p[i]*p[i]<=n;++i) {
            ans-=self(self,n/p[i],i);
        }
        return ans;
    };
    auto phi=[&](auto self,int n,int k) -> int {
        if(n<=p[k]) return 0;
        int ans=g1[idx(n)]-g0[idx(n)]-(fp1[k]-k);
        for(int i=k+1;i<=tot&&p[i]*p[i]<=n;++i) {
            for(int c=1,v=p[i];v<=n;++c,v*=p[i]) {
                ans+=(v-v/p[i])*(self(self,n/v,i)+(c>1));
            }
        }
        return ans;
    };
    printf("%lld %lld\n",phi(phi,n,0)+1,mu(mu,n,0)+1);
}
signed main() {
    int T;
    scanf("%lld",&T);
    while(T--) solve();
    return 0;
}

* III. [HDU-6417] Rikka with APSP

Problem Link

題目大意

給定一張 $n$ 個點的無向完全圖，$i,j$ 之間的邊權定義為最小的正整數 $x$ 使得 $ijx$ 是完全平方數。

求 $\sum_{1\le i\le j\le n}\mathrm{dist}(i,j)$ 的值，其中 $\mathrm{dist}(i,j)$ 表示 $(i,j)$ 之間的最短路長度。

資料範圍：$n\le 10^{10}$。

思路分析

先把每個數 $k$ 寫成標準分解形式 $\prod p_i^{c_i}$，注意到我們只關心 $c_i\bmod 2$ 的值，因此可以把所有的這些 $c_i\bmod 2$ 寫成一個二進位制數 $B_k$。

那麼 $w(u,v)$ 就是 $B_u\oplus B_v$ 中所有 $1$ 對應的 $p_i$ 的乘積，而 $\mathrm{dist}(u,v)$ 就是 $B_u\oplus B_v$ 中所有 $1$ 對應的 $p_i$ 的乘積，先逐步去掉 $B_u\cap(B_u\oplus B_v)$ 中的 $1$，再逐步加入 $B_v\cap(B_u\oplus B_v)$ 中的 $1$，容易證明整個過程中的數值 $\le \max(u,v)$。

因此考慮拆貢獻，對於某個質數 $p$，記 $f(p)$ 為 $[1,n]$ 中含有奇數個 $p$ 的數的數量，那麼 $p$ 對答案的貢獻就是 $p\times f(p)\times(n-f(p))$。

容易發現 $f(p)=\left\lfloor\dfrac{n}{p}\right\rfloor-\left\lfloor\dfrac{n}{p^2}\right\rfloor+\left\lfloor\dfrac{n}{p^3}\right\rfloor-\dots$，考慮根號分治，對於 $p\le \sqrt n$ 的 $p$，暴力計算，時間複雜度為 $\mathcal O(\pi(\sqrt n)\log n)=\mathcal O(\dfrac{\sqrt n}{\log n}\log n)=\mathcal O(\sqrt n)$。

對於 $p>\sqrt n$ 的 $p$，$f(p)=\left\lfloor\dfrac np\right\rfloor$，用整除分塊的技巧列舉 $\left\lfloor\dfrac np\right\rfloor=k$，那麼 $f(p)\times (n-f(p))$ 容易計算，那麼我們只要求某個特定區間 $[l,r]$ 之中的質數之和，考慮整除分塊的過程，注意到 $l-1,r$ 都能寫成某個 $\left\lfloor\dfrac nt\right\rfloor$ 的形式。

因此我們只預處理在所有 $1\sim \left\lfloor\dfrac nt\right\rfloor$ 之間質數之和，在整除分塊的過程中就可以 $\mathcal O(1)$ 回答對應的區間素數數量了。

容易發現這是一個標準的 Min-25 篩的第一部分，因此直接用 Min-25 篩處理即可。

這一部分的時間複雜度為 $\mathcal O(\dfrac{n^{3/4}}{\log n})$。

但是我們還漏掉了一些貢獻，對於 $u\times v$ 是完全平方數但 $u\ne v$ 的數對 $(u,v)$ 對答案是有 $1$ 的貢獻的，因此我們要求，有多少對 $(u,v)$ 滿足 $1\le u<v\le n$，且 $u\times v$ 是完全平方數，考慮構造一個等比數列 $\{v,\sqrt{uv},u\}$，其公比 $<1$ 且所有元素都是 $[1,n]$ 之間正整數，我們只需要對這樣的等比數列計數即可。

設公比最簡分數為 $q/p$，列舉分母 $p$，此時分子 $q<p$ 且 $\gcd(p,q)=1$，因此合法的 $q$ 恰好有 $\varphi(p)$ 種選擇。

然後考慮 $v$，注意到 $u=\dfrac {vq^2}{p^2}$，因此 $p^2\mid v$，這樣的 $v$ 共有 $\left\lfloor\dfrac n{p^2}\right\rfloor$ 種，容易證明這樣的每一對 $(p,q,v)$ 都和我們要計數的 $(u,v)$ 構成雙射。

因此這部分的答案就是 $\sum_{p=2}^{\sqrt n} \varphi(p)\times \left\lfloor\dfrac n{p^2}\right\rfloor$，預處理出 $\varphi(1)\sim \varphi(\sqrt n)$ 即可做到 $\mathcal O(\sqrt n)$。

將答案分成三部分分別處理即可。

第三部分的另一種處理方法

對於這樣乘積是完全平方數的數對計數，最樸素的思路就是列舉 $u,v$ 在各自除掉最大的完全平方因子後得到的數 $k$，列舉這樣的 $k$，$u,v$ 的選擇方案數就是 $g(\lfloor n/k\rfloor)=\binom{\sqrt{\lfloor n/k\rfloor}}2$。

考慮刻畫 $k$ 以寫出求和式，顯然 $k$ 中無平方因子，則 $\mu(k)\ne 0$，因此滿足條件的 $k$ 可以用 $\mu^2(k)=1$ 表示，因此這部分的貢獻就是：

\[\sum_{k=1}^n \mu^2(k)\binom{\sqrt{\lfloor n/k\rfloor}}2 \]

考慮整除分塊列舉 $\lfloor n/k\rfloor$，原問題變成統計 $\mu^2(k)$ 在所有 $\lfloor n/t\rfloor$ 處的字首和，注意到 $\mu^2(k)$ 是積性函數，因此可以用非遞迴版本的 Min-25 篩求出。

由於 $\mu^2(k)$ 函數的特殊性，在遞推 $S(n,k)$ 的時候只需要處理 $c=1$ 的情況，因此該演演算法的執行效率和前一個演演算法的執行效率差別並不大。

兩種做法時間複雜度均為 $\mathcal O(\dfrac{n^{3/4}}{\log n})$。

程式碼呈現

實現方式一：

#include<bits/stdc++.h>
#define int long long
#define LL __int128
using namespace std;
const int MAXN=2e5+1,MOD=998244353;
bool isco[MAXN];
int p[MAXN],tot,val[MAXN],fp[MAXN],phi[MAXN];
int idx1[MAXN],idx2[MAXN],g[MAXN];
inline void solve() {
    int n,m=0;
    scanf("%lld",&n);
    int B=sqrt(n);
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l),val[++m]=n/l;
        if(val[m]<=B) idx1[val[m]]=m;
        else idx2[n/val[m]]=m;
    }
    auto idx=[&](int v) { return (v<=B)?idx1[v]:idx2[n/v]; };
    for(int i=1;i<=m;++i) {
        int k=val[i]%MOD;
        g[i]=(k*(k+1)/2+MOD-1)%MOD;
    }
    for(int k=1;k<=tot;++k) {
        for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {
            g[i]=(g[i]+MOD-p[k]*(g[idx(val[i]/p[k])]+MOD-fp[k-1])%MOD)%MOD;
        }
    }
    int ans=0;
    for(int i=1;i<=tot&&p[i]*p[i]<=n;++i) {
        int s=0;
        for(int v=p[i],c=1;v<=n;v*=p[i],c=-c) s+=c*(n/v);
        ans=(ans+(n-s)%MOD*(s%MOD)%MOD*p[i]%MOD)%MOD;
    }
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l); int k=n/l;
        if((LL)l*l>n) {
            ans=(ans+(n-k)%MOD*(k%MOD)%MOD*(g[idx(r)]+MOD-g[idx(l-1)])%MOD)%MOD;    
        }
    }
    for(int i=2;i*i<=n;++i) ans=(ans+(n/(i*i))%MOD*phi[i]%MOD)%MOD;
    printf("%lld\n",ans);
}
signed main() {
    for(int i=2;i<MAXN;++i) {
        if(!isco[i]) p[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot&&p[j]*i<MAXN;++j) {
            isco[p[j]*i]=true,phi[p[j]*i]=phi[i]*(p[j]-1);
            if(i%p[j]==0) { phi[p[j]*i]=phi[i]*p[j]; break; }
        }
    }
    for(int i=1;i<=tot;++i) fp[i]=(fp[i-1]+p[i])%MOD;
    int T;
    scanf("%lld",&T);
    while(T--) solve();
    return 0;
}

實現方式二：

#include<bits/stdc++.h>
#define int long long
#define LL __int128
using namespace std;
const int MAXN=2e5+1,MOD=998244353;
bool isco[MAXN];
int p[MAXN],tot,val[MAXN],fp1[MAXN];
int idx1[MAXN],idx2[MAXN],g0[MAXN],g1[MAXN],S[MAXN];
inline void solve() {
    int n,m=0;
    scanf("%lld",&n);
    int B=sqrt(n); tot=0;
    for(int i=2;i<=B;++i) {
        if(!isco[i]) p[++tot]=i;
        for(int j=1;j<=tot&&p[j]*i<=B;++j) {
            isco[p[j]*i]=true;
            if(i%p[j]==0) break;
        }
    }
    for(int i=1;i<=tot;++i) fp1[i]=(fp1[i-1]+p[i])%MOD;
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l),val[++m]=n/l;
        if(val[m]<=B) idx1[val[m]]=m;
        else idx2[n/val[m]]=m;
        int k=val[m]%MOD;
        g0[m]=k-1,g1[m]=(k*(k+1)/2+MOD-1)%MOD;
    }
    auto idx=[&](int v) { return (v<=B)?idx1[v]:idx2[n/v]; };
    for(int k=1;k<=tot;++k) {
        for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {
            g0[i]=(g0[i]+MOD-(g0[idx(val[i]/p[k])]-(k-1)))%MOD;
            g1[i]=(g1[i]+MOD-p[k]*(g1[idx(val[i]/p[k])]+MOD-fp1[k-1])%MOD)%MOD;
        }
    }
    for(int i=1;i<=m;++i) S[i]=g0[i];
    for(int k=tot;k>=1;--k) {
        for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {
            S[i]=(S[i]+MOD+S[idx(val[i]/p[k])]-k)%MOD;
        }
    }
    for(int i=1;i<=m;++i) S[i]=(S[i]+1)%MOD;
    int ans=0;
    for(int i=1;i<=tot&&p[i]*p[i]<=n;++i) {
        int s=0;
        for(int v=p[i],c=1;v<=n;v*=p[i],c=-c) s+=c*(n/v);
        ans=(ans+(n-s)%MOD*(s%MOD)%MOD*p[i]%MOD)%MOD;
    }
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l); int k=n/l;
        if((LL)l*l>n) {
            ans=(ans+(n-k)%MOD*(k%MOD)%MOD*(g1[idx(r)]+MOD-g1[idx(l-1)])%MOD)%MOD;    
        }
    }
    for(int l=1,Q=B+5,r;l<=n;l=r+1) {
        r=n/(n/l); int k=n/l;
        while(Q*Q>k) --Q; //Q = sqrt(k)
        ans=(ans+(Q*(Q-1)/2)%MOD*(S[idx(r)]+MOD-S[idx(l-1)])%MOD)%MOD;
    }
    printf("%lld\n",ans);
}
signed main() {
    int T;
    scanf("%lld",&T);
    while(T--) solve();
    return 0;
}